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ABSTRACT 


/ Autocorrelation  Pole-Zero  modeling  identifies  the  parameters  of 
a rational  transfer  function  H(z)  whose  short  time-lag 
autocorrelations  either  exactly  match  (Autocorrelation  partial 
Realization)  or  closely  approximate (Autocorrelation  Prediction^  those 
of  a given  spectrum.  As  a result,  the  spectrum  of  the  H(z)  obtained 
from  either  method  approximates  the  gross  structure  of  the  given 
spectrum.  Autocorrelation  Partial  Realization  (APR)  uses  the 
Pade  approximation  to  determine  the  denominator  coefficients  of  H(z). 
To  compute  the  numerator  coefficients  of  H(z),  APR  uses  an 
iterative  algorithm  such  as  Fejer's  or  Newton-Raphson 's . \In 
contrast,  Autocorrelation  Prediction  (AP)  uses  only  Linear 
Prediction  (LP)  to  determine  both  the  denominator  and  numerator 
coefficients.  Therefore,  once  the  autocorrelation  function  of  the 
given  spectrum  is  known,  AP  uses  only  linear  operations  and 
no  Fourier  Transformations  to  determine  the  parameters  of  H(z). 
Moreover,  the  resulting  rational  transfer  function  is  guaranteed  to 
be  minimum  phase  and  consequently  stable.  AP  can  also  automatically 
determine  the  least  (parsimonious)  denominator  and  numerator  orders 
required  to  model  efficiently  a given  spectral  envelope. 

t-  A dynamic  filtering  process,  based  on  Wiener  filtering  and 
Autocorrelation  Prediction,  was  developed  to  suppress  the  background 
noise  from  degraded  speech.  More  important,  using  AP,  a Linear 


Predictive  Vocoder  was  integrated  into  the  so  called  "Pole-Zero 
Vocoder" (PZV) . Computer  simulations  of  both,  the  dynamic  filtering 
process  and  the  PZV  were  successfully  used  in  speech  processing. 
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CHAPTER  1 


INTRODUCTION 


1 . 1 Problem  Presentation 

Spectral  Pole-Zero  modeling  has  been  the  focus  of  many  research 
efforts  in  recent  years.  The  success  of  Linear  Prediction  [27,29]  in 
All-Pole  modeling  of  the  spectral  envelope  [21,23]  has  encouraged  the 
search  for  a technique  of  comparable  success  in  spectral  envelope 
Pole-Zero  modeling.  On  the  other  hand,  a high  order  All-Pole  model 
is  required  to  model  the  spectral  envelope  having  deep  valleys.  This 
is  so  because  a large  number  of  poles  are  required  to  approximate  a 
small  number  of  zeros  represented  by  these  valleys.  Therefore,  to 
model  a spectral  envelope  having  deep  valleys,  the  Pole-Zero  model 
requires  fewer  parameters  than  the  All-Pole  model.  This  advantage  of 
the  Pole-Zero  model  over  the  All-Pole  model  is  appreciated  in  data 
compression . 

Different  methods  of  Pole-Zero  modeling  have  been  proposed  for 
matching  the  envelope  or  some  smoothed  version  of  a given  spectrum. 
These  methods  may  be  divided  into  two  groups.  The  first  group 
encompasses  those  which  estimate  the  pole  and  zero  parameters 
simultaneously . In  contrast,  the  methods  encompassed  in  the  second 
group  determine  the  pole  and  zero  parameters  separately . The  methods 
in  the  first  group  commonly  face  the  problem  of  solving  a system  of 
non-linear  equations  in  terms  of  the  sought-for  parameters.  Whereas, 
those  in  the  second  group  commonly  use  variations  of  Linear 
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Prediction  to  determine  the  pole  and/or  zero  parameters.  Linear 
Prediction  itself  leads  to  a system  of  linear  equations  solvable  by 
fast  recursive  algorithms  [ 27 ,29 .Appendix  A]. 

Some  recent  f reouencv  domain  methods  of  the  second  group  are 
mentioned  in  the  following.  Cepstral  Prediction  [6, 33.40]  and 
Homomorphic  Prediction  [20,21]  both  use  Homomorphic  deconvolution 
[31]  to  initially  smooth  the  given  spectrum  . Then  each  of  them  uses 
its  own  approach  to  model  the  cepstrally  smoothed  spectrum  with  a 

Pole-Zero  model.  Cepstral  Prediction  is  successful  when  the 
resulting  smoothed  spectrum  is  some  rational  spectrum.  Homomorphic 

Prediction  was  reported  successful  when  applied  to  the  short-time 
spectrum  of  natural  speech  [20,21].  Makhoul  proposed  a method  [27] 
in  which  a variation  of  Linear  Prediction  [10,41]  is  used  to 
determine  the  pole  parameters.  The  zero  parameters  are  then 
estimated  by  Inverse  Linear  Prediction  [27].  This  procedure 
basically  applies  Linear  Prediction  to  reciprocal  of  the  smoothed 
ratio  of  the  given  spectrum  to  the  computed  All-Pole  spectrum. 

Autocorrelation  Partial  Realization  and  Autocorrelation 
Prediction,  described  in  chapter  3,  are  new  techniques;  modeling  the 
envelope  of  a given  spectrum  with  a Pole-Zero  model.  Both  techniques 
determine  the  pole  and  zero  parameters  separately . Also,  both  are 
autocorrelation  domain  methods.  In  other  words,  once  the 
autocorrelation  function  of  the  given  spectrum  is  known,  all  the 
operations  performed  are  in  the  autocorrelation  domain  and  no  Fourier 
Transformation  is  required. 

Autocorrelation  Prediction  was  successfully  used  for  Pole-Zero 
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modeling  in  two  applications:  Wiener  filter  spectral  matching  and 
natural  speech  short-time  spectral  matching. 

1 . 2 Chapter  Summaries 

The  main  body  of  the  dissertation  is  devoted  to  development  of 
Pole-Zero  modeling  techniques  and  their  applications.  The  first  part 
of  the  main  body  develops  new  techniques  for  spectral  Pole-Zero 
modeling  while  the  last  part  applies  these  techniques  to  natural 
speech  processing.  To  make  the  dissertation  more  self-contained,  a 
collection  of  Appendices  is  also  added.  These  Appendices  provide  the 
basic  algorithms  and  concepts  used  in  developing  the  above 
techniques . 

Chapters  2 and  3 comprise  the  first  part,  with  chapter  2 
providing  the  background  needed  for  chapter  3-  The  new  Pole-Zero 
modeling  techniques,  APR  and  AP,  are  derived  in  chapter  3-  Chapters 
4 and  5 comprise  the  last  part,  presenting  the  applications  of  these 
techniques.  A dynamic  filtering  process  that  suppresses  the 
background  noise  from  degraded  natural  speech  is  described  in  chapter 
4.  Chapter  5 presents  the  "Pole-Zero  Vocoder"  (PZV);  an 
analysis-synthesis  process  of  natural  speech  based  on  Autocorrelation 
Prediction . 

The  author's  majors  contribution  are:  the  Pole-Zero  modeling 
techniques  derived  in  chapter  3.  the  improved  dynamic  filtering 
process  described  in  chapter  4,  and  the  integration  of  the  Linear 
Prediction  Vocoder  into  the  Pole-Zero  Vocoder  presented  in  chapter  5. 
Introducing  a new  spectral  flatness  measure,  section  2.2,  and 
preliminary  work  on  a modified  Pole-Zero  Vocoder  to  account  for 
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\ 


background  noise  in  degraded  natural  speech,  section  5.5,  are  some  or 
the  author's  minor  contributions. 


CHAPTER  2 


PARAMETRIC  MODELING  OF  SPECTRA 


2. 1 Introduction 

The  identification  of  a parametric  model  whose  spectrum 
approximates  a given  spectrum  by  minimizing  some  distance  measure  is 
the  general  theme  of  the  chapter.  A new  distance  measure  having  an 
upper  bound  proportional  to  a well-known  distance  measure  E is 
defined.  Two  types  of  parametric  models,  namely  All-Pole  and 
Pole-Zero,  are  focused  upon.  To  show  the  effect  of  the  order  of  the 
All-Pole  model  on  the  approximation,  the  relationship  between  the 
autocorrelation  functions  of  the  model  and  the  given  spectra  is 
derived.  Also,  the  dependence  of  the  minimum  of  the  distance  measure 
E on  the  order  of  the  All-Pole  model  is  analyzed.  Finally,  it  is 
shown  that  minimizing  E to  estimate  the  parameters  of  the  Pole-Zero 
model  leads  to  a system  of  non-linear  equations,  in  contrast  to  the 
system  of  linear  equations  for  the  All-Pole  model. 

2.2  Spectral  Matching  by.  Inverse  Filtering 

The  spectrum  of  the  scaled  model  G H(z)  can  be  matched  to  a 
given  spectrum  | S (a) ) | 2 by  requiring  that  when  S(z)  is  filtered  by  the 
inverse  model  H *(z),  the  spectrum  of  the  output  E(z)  is  flat . 


Figure  2-1 . 
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S(2) 


E ( 2 ) 


Figure  2-1 . 

Inverse  filtering. 

The  model  scale  factor  squared 

G2 

is  equal  to  the 

average 

of  the 

inverse  filter  output  spectrum 

E, 

i.e. 

7T 

2 A 1 

G = E “ 2? 

| E (a>)  | 2 do). 

(2- 

■1) 

-1 

T 

, . . .2 

2 

where,  from  Figure  2-1,  the  output  spectrum  |E(u))|  is  equal  to 


| E(uj)  | 2 « | S(oj)  | 2 | H-1  (oj)  | 2. 


(2-2) 


A distance  measure  for  the  flatness  of  the  output  spectrum  |E(w)|2 
can  be  quantified  as:  The  average  distance  of  the  output  spectrum 

2 

|e(oj)|  from  its  average  E,  i.e. 

it. 


| E(u))  | 2 


E 


dun 


(2-3a) 


-IT 


The  smaller  the  distance  measure  e,  the  flatter  the  output  spectrum 
2 

|e(uj)|  , and  consequently  there  is  a closer  match  between  the  scaled 

2 2 
model  spectrum,  |g  H(w)|  , and  the  given  spectrum  |s(u>)|  . This 

non-negative  distance  measure  has  zero  value  when  the  output  spectrum 

is  constant,  that  is  E(u>)sE  for  -tKuKtt. 
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Properties  of  th£  Distance  Measure  r_.  Using  (2-2),  and  (2-1)  in 
(2- 3a),  gives  another  expression  for  e,  i.e. 

IT 

e - J-  E I S (w)  I 2 H-1(w)  I 2 - | S-1(co)  | 2 dw.  (2-3b) 

— 7T 

2 

The  relation  (2-3b)  shows  that,  for  a given  spectrum  |S(to)|  , the 
optimum  inverse  model  H *(oj)  minimizes  e in  weighted  least  mean 
amplitude  sense.  The  weighting  function  is  proportional  to  the  given 
spectrum  |S(to)|  . Hence,  from  (2-3b),  the  scaled  inverse  model 
spectrum  H ^(w) | 2 approximates  |S  1(w)|2  more  accurately  at  those 
frequencies  where  the  given  spectrum  |s(uj)|  has  its  peaks  rather 
than  valleys. 

An  upper  bound  for  e may  be  found  by  using  the  fact  that  the 
absolute  value  of  the  difference  has  the  least  upper  bound  equal  to 
the  sum  of  the  absolute  values,  i.e. 

| E (co)  | 2 - E < | E(u>)  | 2 + E.  (2-4) 

From  (2-1)  it  is  clear  that  E is  always  positive  except  for  the 
trivial  case,  E(u)  s 0 for  all  ui  in  which  case  E = 0.  Thus  for 
non-trivial  G (uu ) , the  equality  in  (2-H)  holds  only  at  those 
frequencies  where  E(w)=0.  For  non-trivial  E (oo ) , from  (2-4)  and 
(2-3a),  we  obtain 
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n 

C < k { 1 |E(^)|2  + E]du>.  (2-5a) 

—IT 

Performing  the  integration  in  (2-5a)  and  using  (2-1)  gives  the 
following  upper  bound  for  e. 


e < 2E. 


(2-5b) 


Thus  the  flatness  distance  measure  e has  an  upper  bound  proportional 
to  E,  the  average  of  the  inverse  filter  output  spectrum.  Note  that 
the  model  Hq(z)  which  minimizes  the  upper  bound  2E  does  not 
necessarily  minimize  the  flatness  distance  measure  e . It  is 
guaranteed,  nevertheless,  that  the  flatness  distance  measure  for 

Hq(z)  does  not  exceed  the  smallest  upper  bound  2Emln.  Since 
identification  of  the  model  which  minimizes  e is  complicated,  we  seek 
only  the  optimal  model  which  minimizes  the  upper  bound . 2E,  of  the 

\ 

distance  measure  e. 


Other  Interpretations  For  g.  Using  (2-2)  in  (2-1)  gives  another 
interpretation  for  E,  i.e. 

7T 

E “ iS^)|2/  |H(ou)f  2da).  (2_6) 

-IT 


\ 

\ 

\ 


Thus  E is  the  average  of  the  ratio  of  the  given  spectrum  to  the  model 
spectrum.  Hence,  for  the  optimal  model  H(z),  the  average  of  the 
ratio  of  the  given  spectrum  to  the  model  spectrum  is  minimum. 

According  to  Parseval's  Theorem  [32,36],  E given  by  (2-1),  is 
also  the  energy  in  the  inverse  filter  output  signal  E(z). 
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2 . 3 All-Pole  Spectral  Matching 

Let  H(z)  in  Figure  (2-1)  of  section  2.2  be  a stable  All-Pole 
model  defined  as 


H(z) 


1 

A(z) 


1 


l a (k)z 
k=0 


a(0)  - 1.  (2-7) 


1 

To  match  the  spectrum  of  the  scaled  H(z)  to  the  given  spectrum  |s(oo)|  , 
the  prediction  coefficients  (a(k)}  are  computed  to  flatten  the  output 
spectrum  | E (to ) | 2 by  minimizing  the  output  energy  E,  [ 2U , 25 ] . The 

spectrum  of  the  resulting  scaled  H(z)  is  called  the  Linear 
Prediction  (LP)  Spectrum. 

Using  (2-7)  in  (2-6)  gives 

7T 

E * 2tT  j |A(w)|  do).  (2-8) 

-tr 


E is  minimized  by  setting 


9E 

3a(i) 


0, 


1 < i < M. 


(2-9) 


On  the  other  hand,  the  following  relation  holds 

aiW  ' *7i)  lA<‘)k‘‘W-jk“>] 


M 


- [e'31“(I  .(k).!*-)  + .^{J  .(k).-Jk“)l, 

k-0  k-0 


1 < i < M, 


or 
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a 

3a(i) 


|A(u) 


2 


M 

2^  a(k)Cos(i-k)a), 

k-0 


1 < i < M. 


(2-10) 


From  (2-8),  (2-9)  and  (2-10)  and  interchanging  of  summation  and 
integration,  one  obtain 


M i f 

Loa<k)  sf 


|S(uj)|  Cos(i-k)aj  doj  = 0, 


1 < i < M. 


(2-11) 


— TT 


We  know  that  the  autocorrelation  function  R(k)  is  the  Inverse  Fourier 
Transform  of  the  spectrum  |S(u))  | 2,  i .e . 


R(k)  - 


|S(w)|^e 


— TT 


(2-12a) 


or  in  a more  simplified  form 
TT 

R(k)  “ JS'j  |s(w)  |2Cos(kw)dw.  (2-12b) 

-TT 

Substituting  ( 2— 1 2b ) in  (2-11)  results  in  the  well-known 
autocorrelation  normal  equations 

M 

l a (k)R(i-k)  - 0,  1 < i < M.  (2-13) 

k-0  ~ ~ 


The  optimum  predictor  coefficients  { a ( k ) } are  obtained  by 
solving  the  system  of  linear  equations  (2-13)  using  Levinson's  [37] 
or  Durbin's  [27]  recursive  algorithms  [Appendix  A].  Using  (2-13), 
(2-12)  in  (2-8)  and  some  simplification  gives  the  minimum  output 
energy  EM 


^ ' G»  = I a(k)R(k). 
k-0 


(2-14) 


Therefore,  from  Section  2.2,  the  spectrum  of  the  scaled  H(z),  i.e. 


|GAH(eJk“) 


I A (ou)  | ' 


I"  .00<rJk“ 

k=0 


2’ 


(2-15) 


is  the  optimal  All-Pole  spectral  match  to  the  given  spectrum  |S(w)|2, 

Figure  2-2.  The  predictor  coefficients  { a ( k ) } in  (2-15)  and  the 

2 

scale  factor  squared  G.  are  computed  from  (2-13)  and  (2-14), 

A 

respectively . 

2.4  Autocorrelation  Function  of  £hg. 

Optimal  All-Pole  Spectrum 

The  relationship  between  the  autocorrelation  functions  of  the 

2 2 o 

optimal  All-Pole  spectrum  Ga/|a(<d)|  and  the  given  spectrum  |S(u>)|zis 
discussed  here.  Rearranging  (2-15)  gives 


ENERGY  (LOG  SCALE) 


FREQUENCY  (HZ) 


(a) 


(b) 


FREQUENCY  (HZ) 
(c) 


Figure  2-2. 

(a) 

(b) 

(c) 


14-Pole  LP  spectra  superimposing  : 

The  given  short-time  spectrum  of  the  voiced  sound  [0]. 
The  given  short-time  spectrum  of  the  unvoiced  sound  [P]. 
The  given  short-time  spectrum  of  the  unvoiced  sound  [f]. 
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|G  H(eJkw)|2 
A 


A(co) 


A (to) 


(2-16) 


Where  indicates  the  complex  conjugate.  Taking  the  Inverse 

Fourier  Transform  of  both  sides  of  (2-16)  results  in  the  difference 
equation  that  the  autocorrelation  function  of  the  scaled  optimal 
All-Pole  model  obeys,  i.e. 

M 

l a (k)R(i-k)  - G>(-i),  (2-17) 

k-0  A 


where 


R(k) 


JL 

2n 


-IT 


iA(a>)  r 


-jkw 

e J dto. 


(2-18a) 


or 


7T  r2 

if  ^A 

R(k)  - -rr  Cos (kto)dto, 

.1  !*<“>! 2 


(2- 18b) 


and 


h(-i) 


A*((o) 


-j  i(o 
e dio. 


(2-18c) 


Since  h(i)  is  casual  and  from  (2-7) 
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h(0)  = 1, 


(2-18d) 


then  for  i>0  (2-17)  reduces  to 

M r 

l a(k)R(i-k)  =0,  i > 0,  (2-19) 

k=0 


l a(k)R(k)  = G2  1=0  (2-20) 
k=0  A 

The  system  of  equations  given  by  (2-19)  for  1 <i<M  and  (2-20)  has  the 
same  form  as  the  system  of  equations  given  by  (2-13)  and  (2-20). 
Therefore,  these  two  systems  of  equations  have  identical  solutions, 
i .e. 


R(| 1 1 ) - R(|i|). 


0 < i < M 


(2-21) 


Thus,  matching  the  spectrum  of  a scaled  All-Pole  model  of  order  M to 
a given  spectrum  is  equivalent  to  finding  an  All-Pole  model  whose 
scaled  autocorrelation  function  exactly  matches  that  of  the  given 
spectrum  for  the  first  M+1  time-lags. 

Equations  (2-21)  show  that  the  first  M+1  autocorrelations  of  the 
matching  All-Pole  spectrum  G^/|A(co)|  are  exactly  the  same  as  those 
of  the  given  spectrum  |S(u>)|  . The  rest  of  the  autocorrelation 
R(k)  for  | k | >M  are  the  extrapolation  of  R(k)  for  1 <k<M , and  are 
determined  recursively  from  (2-19)  using  the  optimum  predictors  {a(k)l 
|s((i))|  and  Ga/|A(u>)  | are  the  Fourier  Transforms  of  the 
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autocorrelation  functions  R(k)  and  R(k),  respectively.  Therefore, 
increasing  the  order  M of  the  matching  All-Pole  spectrum  increases 
the  range  of  time-lags  over  which  (2-21)  holds,  resulting  in  a closer 

o 2 2 

match  of  the  G7/|A(gj)|  to  | S (to ) | . Hence,  as  we  obtain 
A 


|A»(a))|  l a(i)e_jka) 

k=0 


I S (u) , 


(2-22) 


The  G^,/  ( | 2 can  be  considered  the  All-Pole  representation  of  the 

2 

given  spectrum  |s(w)|  [35]. 


2 . 5 Analysis  of  the  Minimum  Energy  E 

M 

The  minimum  energy  EM,  given  by  (2-114),  of  the  inverse  filter 
output  is  a monotonically  decreasing  function  of  the  All-Pole  model 
order  M.  Figure  2-3  shows  the  normalized  minimum  energy 


vM  - Em/R(0), 


(2-23) 


as  a function  of  the  order  M.  It  can  be  shown  that  0<vMO  for  all  M 
[25,42].  The  minimum  energy  decreases  sharply  as  the  order  M is 
incresed  up  to  some  order  Mp.  For  increases  of  the  All-Pole  model 
order  beyond  Mp,  the  minimum  energy  EM  decreases  slightly.  Mp  is 
referred  to  as  the  parsimonious  (most  economical)  order.  By  the  time 
the  order  M has  reached  Mp,  the  scaled  All-Pole  model  spectrum 
matches  all  the  spectral  envelope  peaks  of  the  given  spectrum.  The 
parsimonious  number  of  poles  Mp.  therefore,  depends  on  the  number  of 
spectral  envelope  peaks  of  the  given  spectrum. 


Figure  2-3.  (a)  Normalized  minimum  energy  curves  for  speech 

short-time  spectra  (voiced  and  unvoiced). 

(b)  Akaike's  Information  Criterion  I(M)  for  voiced 
speech.  The  parsimonious  order  Mp  occurs  at  the 
global  minimum  of  I(M),  shown  by  the  arrow  at  M=10. 

(c)  Akaike's  Information  Criterion  I(M)  for  unvoiced 
speech.  The  parsimonious  order  Mp  occurs  at  the 
global  minimum  of  I(M),  shown  by  the  arrow  at  M=3. 
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To  estimate  Mp  for  a given  spectrum,  as  the  order  M of  the 
All-Pole  model  increases,  one  might  monitor  the  normalized  minimum 
energy  given  by  (2-23).  This  monitoring  can  be  done  automatically 
using  Durbin's  recursive  algorithm,  a simplification  of  Levinson's 
algorithm  [Appendix  A].  Durbin's  recursive  algorithm  computes  the 
minimum  energy  EM  for  successive  orders  as  a by-product  of  the 
calculation  of  the  predictor  coefficients. 

Another  approach  to  estimate  Mp,  when  the  length  of  the  time 
sequence  that  generates  the  given  spectrum  is  known,  was  proposed  by 
Akaike  [ 1 ,2,3,4, 15] . Makhoul  applied  this  approach  to  estimate  Mp  in 
speech  spectral  matching  [26].  In  this  approach,  the  so-called 
Akaike's  Information  Criterion  I(M),  given  by  (2-24)  below,  has  its 
global  minimum  at  Mp. 

2M 

I (M>  = L°8  vm  + ^ • (2-24) 

The  W in  (2-24)  is  the  time  domain  window  length  and  the  constant  c, 
0<c<1,  accounts  for  the  effective  length  of  the  time  domain  window 
(Makhoul  reports  c = o.4  for  a Hamming  window),  Figure  2-3. 

To  match  the  spectral  envelope  peaks  and  valleys  of  the  given 
spectrum,  a large  increase  in  the  order  M beyond  Mp  is  required. 
This  is  so,  because  the  spectral  zeros  represented  by  the  deep 
valleys  require  a large  number  of  poles  to  approximate  them.  To 
avoid  the  large  increase  in  M beyond  Mp,  zeros  are  introduced 
explicitly  into  the  model.  In  other  words,  the  given  spectrum  is 
matched  with  the  spectrum  of  a Pole-Zero  model  rather  than  an 


All-Pole  model. 
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2-6  Pole-Zero  Spectral  Matching  [ 27 1 

Let  H(z)  in  figure  (2-1)  of  section  2.2  be  a Pole-Zero  model 


defined  as 


H(z) 


= B(z) 


A(z) 


L 

l b(i)z 
1=0 
M 

I a(i)z' 

i=0 


(2-25) 


To  match  the  spectrum  of  scaled  H(z)  to  the  given  spectrum  ls(to)|2, 
the  optimal  pole  predictors  (a(k)}  and  the  zero  predictors  (b(k)}  are 
computed  by  minimizing  the  output  energy  E,  [27],  Using  (2-25)  in 
(2-6)  gives 


E = 


2 lA(u))|2 
! B Cco)  | 2 


dw. 


(2-26) 


E is  minimized  by  setting 
9E 

3a(i)  ’ 1 1 1 1 M,  (2-27a) 


9b(i)  °’  1 1 i 1 L.  (2-27b) 

On  the  other  hand,  similar  to  (2-10),  the  following  relation  holds: 


9 

3b(i) 


I B (to)  |2 


L 

2 I b (k)Cos (i-k)w, 
k=0 


1 < i < L. 


(2-28) 


We  also  define 
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Vk) 


7T 


2tt 

-7T 


I S(uj)  | 2 


^ Cos(kw)dw, 

|B(co)|26 


(2-29) 


where  a and  p are  positive  integers.  Using  (2-26),  (2-10)  and  (2-2 9) 
in  (2-27a)  and  some  simplification  results  in 

M 

1 a(k)R  (i-k)  =0,  1 < i < k.  (2-30a) 

k=0 

Similarly,  using  (2-26),  (2-28)  and  (2-29)  in  (2-27b)  and  some 

simplifying  results  in 

L 

l b(k)R  (i-k)  =0,  1 < i < L.  (2-30b) 

From  (2-29)  it  is  clear  that  R^(k)  is  a function  of  the  {b(i)J 
whereas  i3  a function  of  both  {a(k)}  and  (b(k)}. 

Consequently,  the  system  consisting  of  equations  (2-30a)  and  (2-30b) 
is  non-linear  in  terms  of  { a ( k ) } and  (b(k)}.  To  solve  this  system  of 
non-linear  equations . an  iterative  scheme  can  be  used  [27,36]-  An 
iterative  scheme,  however,  brings  about  its  own  problems,  such  as 
convergence,  stability,  and  high  rates  of  computation. 

To  avoid  these  problems,  in  the  next  chapter  we  seek  suboptimal 
Pole-Zero  models  whose  parameters  are  partially  or  totally  solutions 
to  systems  of  linear  equations. 


CHAPTER  3 


\ 


AUTOCORRELATION  FUNCTION  POLE-ZERO  MODELING 
3 • 1 Introduction 

Identif ication  of  a Pole-Zero  model  whose  autocorrelation 
function  approximates  that  of  a given  spectrum  is  the  main  issue  of 
this  chapter.  To  see  the  relationship  between  the  parameters  and  the 
autocorrelation  function  of  a Pole-Zero  model,  the  difference 

equation  governing  the  autocorrelation  function  of  the  Pole-Zero 
model  is  derived.  Two  different  techniques,  "Autocorrelation  Partial 
Realization"  (APR)  and  "Autocorrelation  Prediction"  (Af1),  are 
developed  to  estimate  the  parameters  of  a Pole-Zero  model  whose 
autocorrelation  function  approximates  that  of  a given  spectrum.  In 
either  technique,  the  pole  parameters  and  the  zero  parameters  are 
estimated  separately  . 

APR  uses  the  Pade ' approximation  to  estimate  the  pole  parameters 
and  an  iterative  method  to  estimate  the  zero  parameters  by  solving  a 
system  of  non-linear  equations.  In  contrast,  AP  uses  Linear 
Prediction  to  estimate  both  the  pole  and  the  zero  parameters.  The 
spectral  interpretation  of  AP  is  given  and  the  selection  of  orders  of 
the  Pole-Zero  model  are  discussed.  Finally,  APR  and  AP  are  compared. 

3.2  Autocorrelation  Function  of  Pole-Zero  Model 

To  understand  the  relation  between  the  autocorrelations  of  a 
stable  Pole-Zero  model,  the  difference  equation  that  governs  the 
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Pole-Zero  model  is  derived. 

Consider  the  Pole-zero  model: 


aQ  * 1,  (3-la) 


whose  power  series  representation  is: 


H(z) 


B(z) 

A(z) 


l b.z"1 

i-0 

M -i 

l ajz 

i»0 


H(z)  = l h(k)z_k , 
k«0 


M > r, 

r < 1, 


(3-lb) 


where,  from  (3-1a)  and  ( 3— 1b),  the  coefficients  h(k)  are  obtained 
from  the  following  recursive  formula: 


L M 

h(k)  = l b(i)6(k-i)  - l a(i)h(k-i)  . (3-2) 

i=0  i-1 


The  ( a^  > and  {b^}  are  referred  to  as  pole  predictors  and  zero 
predictors,  respectively.  Multiplying  both  sides  of  (3-1a)  by 
H(1/z)A(z)  gives: 


[H(z)H(l/z) ]A(z)  - B(z)H(l/z).  (3-3) 

Taking  the  Inverse  Z-transform  of  both  sides  of  (3-3)  and  using 
(3-1a)  and  (3-1b)  results  in  the  difference  equation: 

M A L 

l a(i)  R(k-i)  - [ b(i)  h(i-k), 

i-0  i-0 


(3-4a) 
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or 

- L M 

R(k)  - l b(i)h(i-k)  - l a(i)R(k-i)  , (3-4b) 

i-0  i-1 

where  the  model  autocorrelation  function  R(k)  is  defined  by: 

+°°  . 

I R(k)z~k  - H(z)H(l/z)  , r < |z|  < 1/r  , (3-5) 

k— 00 

and  from  ( 3— 1 1> ) the  power  series  representation  of  H(1/z)  is: 

oo  Q 

H(l/z)  = l h(k)zk  = l h(-k)z_k  , |z|  <-.  (3-6) 

k=0  k— °°  r 

Since  h(k)  is  causal,  then  the  difference  equation  (3-4)  is  further 
simplified  for  k > L.  That  is, 

M 

l a(i)R(k-i)  = 0 , k > L . (3-7) 

i-0 

Comparing  (3-4b)  and  (3-2)  reveals  that  for  k > L the  same  pole 

A 

predictors  { } predict  both  R(k)  and  h(k)  from  their  corresponding 
last  values. 

3.3  Aatas.Qrre),3ti9n  Partial  Realization 

The  idea  is  to  find  the  Pole-Zero  model  H(z),  of  the  form 
( 3— la),  whose  autocorrelation  function  Rk  exactly  matches  the 
symmetric  autocorrelation  function  R^ , of  a given  spectrum 
S(z)S(1/z),  for  the  first  N time-lags.  That  is, 
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R|k[  = R|k|  , k = 0,1,2, ...N-l  , (3-8) 

\ 

The  autocorrelation  function  R,  is  the  Inverse  Z-transform  of  the 

k 

given  spectrum  S(z)S(1/z),  i.e. 


I = S(z)S(l/z)  , 

k=-oo 


< 1 . 


(3-9) 


and  similarly,  Rk  is  the  Inverse  Z-transform  of  the  model  spectrum 
H(z)mi/z): 


£ R,z  •’  H(z)H(l/z) 
k=-o° 


B(z)B(l/z) 

A(z)A(l/z) 


r2  < 1- 


(3-10) 


The  idea  is  realized  in  three  steps,  Figure  3-1: 

i.  We  show  that  using  the  Pade'  approximation  [Appendix  B]  on  the 
right  half  of  the  autocorrelation  function  R^  leads  to  a rational 
function  C(z)/D(z)  whose  denominator  is  equal  to  the  sought  for 
denominator  A(z).  Furthermore,  P(z),  the  numerator  of  the  two-sided 
rational  function  C(z)/D(z)+C(1/z)/D(1/z)  is  shown  to  be  the  spectrum 
of  the  desired  numerator  B(z). 

ii.  Some  properties  of  the  polynomial  P(z)  are  discussed,  and  a 
direct  method  for  computing  its  coefficients  is  derived. 

iii.  Some  iterative  algorithm  such  as  Fejer's  [37]  or 

Newton-Raphson 's  [43,  Appendix  D]  is  used  'ro  decompose  the  P(z)  into 
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B(z)B( 1/z) . 

Finally  some  issues  invloved  in  Autocorrelation  Partial 
Realiztion  are  discussed. 


i.  Pade  Approximation . Consider  the  one-sided  sequence  xk 
defined  as: 


*k  = 


k = 0 
k > 0 . 


(3-11) 


Then  the  two-side  power  series  j Rkz_*C  can  be  decomposed  into  the 

k=-°° 

sum  of  the  following  two,  one-sided  power  series: 


-Hr 


l R^z'W  = l x zk  + l x^ 


-k 


k— 00 


k-0 


k=0 


r,  < z < - , 
1 T1 


rl  < i. 


(3-12) 


Furthermore,  the  right  most  power  series  in  (3-12)  is  denoted  by  X(z). 
That  is, 


X(z)  - l 
k-0 


z > r 


1’ 


1 


(3-13) 


Now  , using  the  Pade'  approximation  [Appendix  B],  the  one-sided 
power  series  X(z)  is  approximated  with  the  stable  rational  function 
X(z) : 
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l C4Z  1 

X(Z)  £ Clz)  , k^o  . r xz-k 

U)  D(z)  M , ,LV  ’ 


l dLz 

i-0 


-i  k-0 


dO  " x- 


\z\  > r2> 


r2  " 


(3-14) 


where,  from  (3-14),  the  coefficients  of  the  C(z)/D(z)  power  series 
representation  are  given  by 


*k 


M 

1 


i-0 


c16(k-i) 


diVi  • 


(3-15) 


The  Berlekamp-Massey  [7,30]  or  Trench  [39,41]  recursive  algorithms 
[Appendix  C]  can  be  used  to  perform  the  Pade  approximation,  leading 
to  the  rational  function  X(z)=C(z)/D(z) . As  shown  in  Appendix  C,  the 
Pade^  approximation  X(z)  for  the  power  series  X(z)  has  the  following 
property: 

v - x,  , k - 0,  1,  2,  ...  N-l  , (3-16) 


where 


N - 2M  + 1 . 


(3-17) 


Also  from  (3-17)  and  (3-15)  one  concludes  that 


\ m - l diVi  . 

K * 


k > N . 


(3-18) 
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A 

Therefore,  X(z)=C(z)/D(z)  approximates  the  given  power  series  X(z), 
satisfying  the  equations  (3-16)  and  (3-18).  In  other  words,  the 

A 

first  N coefficients  of  the  power  series  X(z)  are  equal  to  those  of 
the  power  series  X(z)  and  the  rest  are  recursive  extrapolations  of 
these  first  N coefficients.  Thus, 


~D(f)  * lzl  > **x{tvt2)  , 
Max(rltr2)  < 1 . 


We  now  find  a two-sided  rational  spectrum  whose  Inverse 
Z-transform  satisfies  the  equation  (3-8).  Using  (3-19)  in  (3-12)  and 
some  simplifying  operations  results  in 


D(z)D(l/z)  ’ 


Max(r j,r2>  < |z|  < l/Maxtr^.r^)  , 
Max(r1,r2>  < 1 , 

or 


I v 

k— ® 


l V 


D(z)D(l/z) 


D(z)D(l/z)  ’ 


Max(r1,r2)  < |z|  < l/MaxO^.rp  , 
Max(r1>r2>  < 1 , (3-20b) 


where 
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R 


M 


A 

s 


k = 0 

Ho. 


(3-21) 


The  two-sided  rational  spectrum  P(z)/D(z)D( 1/z) , given  in  (3-20b), 
has  the  Inverse  Z-transform  R^  for  which,  from  (3-21),  (3-16)  and 
(3-11),  the  relation  (3-8)  holds.  Since  P(z)/D(z)D( 1/z)  is  the 
spectrum  of  the  desired  Pole-Zero  model  H(z),  then: 


B(z)B(l/z)  _ P(z)  . 
A(z)A(l/z)  D(z)D(l/z) 


(3-22) 


From  (3-22)  the  model  denominator  A(z)  is  equal  to 
A(z)  = D(z)  , (3-23) 

while,  from  (3-22)  and  (3-20)  the  spectrum  of  the  model  numerator 
B(z)  is 

B(z)B(l/z)  - P(z)  - C(l/ z)D(z)  + C(z)D(l/z)  • (3-24) 

ii.  Properties  of  the  P(z) . Before  attempting  to  decompose  the 
polynomial  P(z),  we  examine  some  properties  of  P(z)  and  give  a direct 
method  for  computing  P(z).  From  (3-24)  and  (3-14),  it  can  be  shown 
that  P(z)  is  a symmetric  two-sided  polynomial  with  real  coefficients. 


That  is, 
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P(z)  = l 


P|ulZ 


k r * 


(3-25) 


where 


n- l\  n-K 

J0  Cidi+k  + 1^Q  diCi+k  ’ k = 0,  1,  ...M 


(3-26) 


Therefore  P(z)  has  real  values  on  the  unit  circle  in  the  z plane. 

To  compute  the  polynomial  P(z)  directly,  consider  the  symmetric 
two-sided  power  series  U(z)  defined  as 


+“  , , 
D(z)D(l/z)  l R^z  = l U^z 

l(a-oo  Jq*— oo 


(3-27) 


r1  < |z|  < l/r1  , 


rl  < X- 


From  (3-20b)  we  also  have 


P(z)  - D(z)D(l/z)  l R z~k  , 

k»-oo 


r2  < I Z | < l/r2  , 

r2  < i. 


(3-28) 


Similar  to  P(z),  D(z)D(1/z)  is  also  a symmetric  two-sided  polynomial 
of  order  M.  From  this  property  of  D(z)D(1/z)  and  using  the  equations 
(3-8)  and  (3-17)  in  comparing  (3-27)  with  (3-28),  one  concludes  that 
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k = 0,  1,  . . . M • (3-29) 


(3-29)  shows  P(z),  given  by  (3-25),  is  obtained  by  truncating  the 
power  series  U(z),  i.e. 


P(z) 


M 


I 

k=-M 


M 


l 

k=-M 


(3-30) 


Thus,  the  polynomial  coefficients  pj^j  can  be  computed  directly  by 
convolving  the  autocorrelation  function  of  the  A(z)  coefficients  with 
that  of  the  given  spectrum. 

iii.  Decomposition  of  the  P(z) . To  decompose  P(z)  into 
B(z)B(1/z),  Fejer's  [37]  or  Newton-Raphson ' s iterative  algorithm 
[Appendix  D]  can  be  used.  Fejer's  algorithm  finds  the  2M  roots  of 
the  symmetric  polynomial  P(z)  and  properly  chooses  M of  the  them  to 
construct  the  polynomial  B(z).  The  latter  algorithm  finds  the 
minimum  phase  polynomial  B(z)  whose  spectrum  approximates  P(z)  with 
desired  accuracy. 

Discussion . A few  issues  concerning  Autocorrelation  Partial 
Realization  deserve  attention.  These  are:  a.  the  stability  of  the 
Pole-Zero  model  B(z)/A(z),  b.  the  condition  that  P(z)  should  meet  to 
make  its  decomposition  into  B(z)B(1/z)  possible,  and  finally,  c.  the 
non-linear  methods  required  to  perform  the  decomposition  in  b. 

a.  To  have  an  stable  Pole-Zero  model  B(z)/A(z),  the  rational 
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function  C(z)/D(z)  resulting  from  the  Pade  approximation  should  be 
stable  . We  used  this  assumption  to  show  that  equation  (3-8)  holds 
for  the  Inverse  Z-transform  of  P(z )/D(z )D( 1 /z ) . For  a given  N,  the 
rational  function  C(z)/D(z),  however,  is  not  guaranteed  to  be  stable. 
Because  of  the  equations  (3-11)  and  (3-8),  one  hopes  that  by 
increasing  N and  consequently  M,  one  finally  finds  a stable  C(z)/D(z) 
and  as  a result  a stable  Pole-Zero  model  B(z)/A(z). 

b.  To  decompose  the  P(z)  into  B(z)B(1/z),  the  polynomial  P(z) 
should  be  non-negative  on  the  unit  circle.  Why  this  condition  should 
be  met  is  verified  easily  from  (3-24)  where  B(z)B(1/z)  is 

non-negative  on  the  unit  circle.  On  the  other  hand,  P(z)  is  equal  to 
the  truncated  U(z).  Though  U(z)  is  non-negative  on  the  unit  circle, 
this  is  not  necessarily  true  for  the  truncated  U(z).  Because  of 
(3-27)  and  (3-30),  the  P(z)  can  be  made  non-negative  on  the  unit 
circle,  however,  by  choosing  some  high  order  M. 

c.  Due  to  the  non-linear  equation  (3-24),  the  composition  in 
(b)  requires  an  iterative  method.  Consequently,  the  problems  of 
non-negativeness , convergence  and  the  high  rates  of  computation 
should  be  dealt  with. 


3 . 4 Autocorrelation  Prediction 

The  idea  is  to  find  a minimum  phase  Pole-Zero  model  of  the  form 
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ao  bo  = 1 ’ 

L < M , 


(3-31) 


whose  autocorrelation  function  R(k)  approximates  the  autocorrelation 


function  R(k)  of  a given  spectrum  S(z)S(1/z),  u ing  only  linear 


operations,  i.e. 


R(k)  «=  R(k) , 


(3-32) 


where  R(k)  is  defined  by  (3-9)  and  R(k)  is  defined  by  : 
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(3-33) 


r2  < | z | < l/r2  , 
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Also,  using  (3-31),  the  autocorrelation  functions  Rg(k)  and  R^(k)  in 
(3-33)  are  obtained  from 


rboo 


■ r 

i=0 


bi  bi+|k| 


(3-34) 


M-|k| 

RA(k>  = Lo  ai3i+lkl'  (3"35) 

The  above  idea  is  realized  in  three  steps:  i.  Estimation  of  the 
minimum  phase  polynomial  A(z),  ii.  Computation  of  the  "Residual 


Autocorrelation  Function"  defined  as: 
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where  " ®"  indicates  linear  convolution,  iii.  Estimation  of  the 
minimum  phase  All-Zero  model  B(z)  and  the  gain  G,  Figure  3-2.  At  the 
end,  the  approximation  error,  and  how  Autocorrelation  Prediction 
compares  with  Autocorrelation  Partial  Realization  are  addressed. 

i.  Estimation  of  the  A(z) . The  minimum  phase  polynomial  A(z)  is 

obtained  simply  by  finding  the  All-Pole  model  GA/A(z)  whose 

2 

autocorrelation  function  GA^/A(k)  exactly  matches  that  of  the  given 
spectrum  S(z)S(l/z)  for  the  first  M+1  time-lags,  i.e. 

R(k)  = R1/aOO.  0 < k < M,  (3-37) 


where  the  autocorrelation  function  G.R,  ,.(k)  is  defined  as: 
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(3-38) 
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As  is  shown  in  section  2.3  and  2.4,  the  pole  predictors  { a ^ } are  the 

solution  to  the  system  of  linear  equations  (2-13)  which  also 

2 

guarantees  that  A(z)  is  minimum  phase.  The  gain  squared  GA  is 
obtained  from  (2-14). 


ii.  Residual  Autocorrelation  Function.  After  finding  the  pole 
predictors  {a^},  the  autocorrelation  function  RA(k)  is  computed 
according  to  the  definition  (3-35).  Then,  from  (3-36).  the  Residual 


Figure  3-2.  Autocorrelation  Prediction  block  diagram. 
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Autocorrelation  function  R^(k)  is  obtained  by  simply  convolving  the 

finite  length  autocorrelation  function  R (k)  with  the  autocorrelation 

A 

function  R(k) . 

iii.  Estimation  of  the  B( z ) and  G.  The  estimation  of  the 
All-Zero  model  B(z)  and  the  gain  G is  performed  in  steps  a and  b. 
a.  We  find  a high  order  All-Pole  model  of  the  form 


Q(z) 
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l 

i=0 


qiZ 


-i 


qQ  = 1,  0-39) 

N » L, 


2 

whose  autocorrelation  function  GqR^U)  exactly  matches  the  Residual 
Autocorrelation  function  Rp(k)  for  the  first  N+1  time-lags.  That  is, 


CQ  Rl/q'k>  ‘ Vk)' 


° < f k f < N, 


(3-40) 


2 

where  the  autocorrelation  function  G^R^^Ck) 


is  defined  as 
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(3-41) 


Similar  to  the  estimation  of  the  { a^ } , the  coefficients  { q ^ ) are  the 
solution  to  the  system  of  linear  equations  (2-13)  after  replacing 
{ajJ  by  (q^),  R(k)  for  0_<|k|_<N  by  R^k)  for  0<|k|<N,  and  M by  N,  i.e. 
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N 

l q.R^k-i)  = 0 * 1 < k < N.  (3-42) 


Also,  the  gain  squared  Gg  is  obtained  from  (2-1 4)  after  the  same 
replacement.  That  is, 


cq  ■ L 


(3-43) 


i=0 


Using  (3-39)  in  (3-41)  gives  the  expression  for  the  autocorrelation 
function  Rg(k)  as 


Yk> 
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(3-44) 


After  computing  the  coefficients  {q^}  from  (3-42),  the  equation 
(3-44)  is  used  to  calculate  R^(k). 

b.  Finally,  we  find  the  All-Pole  model  G„/B(z)  whose 

D 

2 

autocorrelation  function  GbR1/B^  exactly  matches  the 
autocorrelation  function  Rg(k)  for  the  first  L+1  time-lags,  i.e. 


gbri/b°°  “ Yk)* 


0 < I k I < L, 


(3-45) 


where  the  G^R^^k)  is  defined  as 
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(3-46) 


Again,  similar  to  the  estimation  of  the  {a^.the  zero  predictors  {Ik} 
are  the  solution  to  the  system  of  linear  equations  (2-13)  after 
replacing  {aj}  by  { bA } , R(k)  for  0<|k|<M  by  Rp(k)  for  0j< | k | _<  1 , and  M 
by  L.  That  is, 


L 

l b R (k-i)  =0,  1 < k < L.  (3-47) 

i=0  y ~ ~ 


Similar  to  (2-13).  the  solution  of  the  system  of  linear  equations 

(3-47)  guarantees  that  B(z)  is  minimum  phase.  Also,  the  gain  squared 
2 

Gg  is  obtained  from  (2-14)  after  the  same  replacement. 
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It  is  shown  in  the  next  section  that  the  gain  G is  equal  to 
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The  identification  process  of  the  Pole-Zero  model  G B(z)/A(z)  is 
summarized  in  Figure  3-2.  Figure  3-3  shows  a speech  short-time 
spectrum  superimposed  by  the  All-Pole  and  Pole-Zero  model  spectra  for 
comparison.  Note  the  advantage  of  the  Pole-Zero  model  over  the 
All-Pole  model  in  matching  the  spectral  envelope  valley  of  the  given 
speech  short-time  spectrum. 

A variation  for  AP  is  obtained  by  computing  the  pole  predictors 
from  a non-symmetric , rather  than  a symmetric,  toepltiz  matrix.  In 
this  way,  the  zeros  in  the  model  are  explicitly  accounted  for.  For 
this  variation,  however,  the  stability  is  no  longer  gauranteed. 

Approximation  Error.  To  show  how  close  the  autocorrelation  of 
the  resulting  Pole-Zero  model  G B(z)/A(z)  approximates  the 
autocorrelation  function  R(k),  the  approximation  error  is  derived. 
From  (3-38),  (3-^1)  and  (3-46),  the  following  similar  equations  are 
obtained 

GA  Rl/A(k)  ® RA(k)  * G^k>.  (3- 50a) 
Gq  Rl/Q(k)  ® Rq(k)  - Gq600*  (3- 50b) 

GB  Rl/B°°  ® RB°°  “ Gb600'  (3- 50c) 

Using  (3-40) , we  define  the  error  autocorrelation  function  ARR(k)  to 
be 
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(b) 

Figure  3-3.  (a)  Speech  short-time  spectrum  superimposed  with 

14-Pole  model  spectrum  (Linear  Prediction). 

(b)  Speech  short-time  spectrum  superimposed  with 

8-Pole  and  6-Zero  model  spectrum  (Autocorrelation 
Prediction). 
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Using  ( 3— 50b ) in  (3-52)  gives 

RQ°°  ® Vk)  = GQ6<k)  + AR^(k)  ® Rq00. 


Similarly,  using  (3-H5),  we  define  the  error  autocorrelation 
ARq(k)  to  be 
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Thus  we  have 
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Using  (3-55)  and  (3-36)  in  (3-53)  and  some  rearranging 
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(3-53) 

function 


(3-54) 
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of  terms 


results  in 
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CB  Rl/B(k)  ® RA(k)  0 R(k)  = GQ6(k)  + [AVk>  ® Rq<k)  - ARQ(k>  ® 

Rj/k)].  (3- 56) 

Finally,  using  (3-50a),  (3-50c)  in  (3-56)  gives 
G0 

R(k)  = -j  RR(k)  ©R1/A(k)  + AR(k),  (3-57) 

GB 

where  the  approximation  error  is 

AR(k)  =-^-RB(k)  ® R1/A(k)  © [ARD(k)  ©RQ(k)  - ARQ(k)  © R^k)]. 

CB 

(3-58) 

For  proper  selection,  discussed  below,  of  the  orders  M,  L,  and  N the 
function  in  the  square  brackets  in  (3-56)  or  (3-58)  has  values 
relatively  close  to  zero  for  the  short  time-lags.  Therefore, 

rn 

R(k)  = RR(k)  © R1/A(k)-  0-59) 

gb 

Taking  Z-transform  of  (3-59)  and  using  (3-9)and  (3-38)  results  in 


l R(k)z“k  = S(z)S(l/z) 
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(3-60) 


Hence,  the  autocorrelation  function  of  the  Pole-Zero  model 
approximates  the  R(k)  at  short  time-lags  and  the  approximation  error 
is  given  by  (3-58).  Note  (3-56)  shows  that  the  autocorrelation 
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function  of  the  inverse  filter  output,  for  the  obtained  Pole-Zero 

2 

model,  differs  from  the  impulse  response  Gfi(k)  by  1 /G  times  the 

B 

’small"  function  in  the  square  brackets 

Comparison  of  APR  and  AP.  There  is  a trade-off  between  the 
Autocorrelation  (AP)  Prediction  and  Autocorrelation  Partial 
Realization! APR)  . In  AP  all  the  operations  are  linear  and  the 
stability  is  guaranteed,  but  the  short  time-lag  autocorrelations  of 
the  given  spectrum  are  approximately  matched . In  contrast,  APR  has 
partially  non-linear  operations  and  the  stability  is  not  guaranteed, 
but  the  short  time-lag  autocorrelations  are  exactly  matched . The  APR 
is  also  theoretically  more  appealing.  AP  and  APR  have,  however,  some 
common  properties.  Once  the  autocorrelation  of  the  given  spectrum  is 
known,  no  Fourier  transformation  is  required  to  estimate  the  model 
parameters , using  either  technique.  More  important,  both  are  well 
suited  for  matching  the  spectral  envelope  of  a given  spectrum  having 
fine  structure,  such  as  speech  short-time  spectrum.  This  is  possible 
because  in  both  techniques  the  model  parameters  are  computed  such 
that  the  short  time-lag  autocorrelations,  representing  the  gross 
structure  of  the  given  spectrum,  are  either  exactly  matched  or 
closely  approximated. 

3.5  Frequency  Interpretation  of  AP 

The  frequency  domain  interpretation  of  the  steps  taken  in 
Autocorrelation  Prediction,  described  in  section  3-4,  and  selection 
of  the  orders  M,  L and  N are  addressed  here,  Figure  3-4. 

The  given  spectrum  S(z)S(1/z)  can  be  thought  of  as  being 
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composed  of  spectral  poles  and  zeros  represented  by  the  spectrum 

peaks  and  deep  valleys.  As  disscussed  in  sections  2.3  asd  2.5, 

2 

fitting  the  All-Pole  speccrum  G /A(z)A(1/z)  to  the  given  spectrum 

A 

S(z)S(1/z),  matches  the  peaks  rather  the  valleys  of  the  spectral 
envelope.  By  selecting  the  order  of  the  All-Pole  model  equal  to  the 
parsimonious  order  M[)f  those  peaks  having  higher  amplitudes  are 
matched,  while  the  deep  valleys  are  left  unmatched.  The  residual 
spectrum  D(z)D(1/z),  which  is  the  given  spectrum  after  removing  the 
estimated  poles, contains  primarily  the  spectral  zeros  represented  by, 
the  deep  valleys  in  the  residual  spectrum.  To  fit  the  All-Zero 
spectrum  B(z)B(1/z)  to  the  residual  spectrum  D(z)D(1/z),  the  All-Pole 
spectrum  G^/B(z)B( 1/z)  is  fitted  to  the  Q(z)Q(1/z);  an  approximation 
for  the  reciprocal  of  the  residual  spectrum  envelope.  To  obtain 


to  the  residual  spectrum  D(z)D(1/z).  The  high  order  N is  required 
because  each  of  the  spectral  zeros  of  the  residual  spectrum  is 
approximated  with  a large  number  of  poles  in  the  model  Gq/Q(z).  The 
lower  bound  for  the  order  N,  therefore,  depends  on  the  number  of 
spectral  zeros  in  the  residual  spectrum  or,  equivalently,  in  the 
given  spectrum  S(z)S(1/z),  and  how  close  these  spectral  zeros  are  to 
the  unit  circle.  Similar  to  the  order  M,  the  order  L may  also  be 
equal  to  the  parsimonious  order  Lp  of  the  All-Pole  spectrum  match  to 
the  spectrum  Q(z)Q(1/z). 

To  reduce  the  computation  of  the  parameter  estimation  when  the 
given  spectral  envelope  has  more  deep  valleys  than  peaks,  the 
Autocorrelation  Prediction  is  applied  to  the  reciprocal  of  the  given 


spectrum  and  then  reciprocal  of  the  resulting  Pole-Zero  model  is 
used.  In  chapter  4 this  approach  is  used  in  Pole-Zero  modeling  of 
Wiener  filter  spectrum. 

3 . b Realization  of  Pole-Zero  Model 

The  Pole-Zero  model  may  be  realized  using  any  of  a variety  of 
methods  [16,17].  We  propose  a method  which  is  more  suitable  for 
realization  of  the  Pole-Zero  model  obtained  by  AP  or  APR. 

Cascade  of  Two  Lattice  Filters . The  Pole-Zero  model  G B(z)/A(z) 
giyen  by  (3-31)  can  be  realized  as  the  cascade  of  the  All-Zero  filter 
B(z)  and  the  All-Pole  filter  1/A(z),  having  the  overall  gain  G. 
Either  of  these  filter  is  implemented  by  a lattice  filter  [19]  which 
is  identified  with  the  so-called  parcor  parameters,  Figure  3-5. 

There  is  a unique  set  of  parcor  parameters  {K^}  for  the  set  of 
predictor  coefficients  [a^]  and  vice  versa  [Appendix  A].  Similarly, 
there  is  unique  set  of  parcor  parameters  { } for  the  zero  predictors 
[bj]  and  vice  versa.  Both  sets  of  parcor  parameters  [Kj } and  {K^ } 
are  obtained  as  a by-product  using  AP,  but  both  should  be  computed 
[Appendix  A ] when  the  Pole-Zero  model  is  identified  using  APR. 


CHAPTER  4 


DYNAMIC  POLE-ZERO  FILTERING 


4 . 1 Introduction 

Quite  often  speech  is  degraded  with  additive  independent 
background  noise,  such  as  interference  of  helicopter  noise  with 
pilot's  speech  and  surface  noise  with  singing  in  the  play  back  of  old 
recordings.  To  enhance  the  speech,  an  input  controlled  time-varying 
filter  is  devised  that  suppresses  the  noise  and  passes  through  the 
speech . 

Speech  and  singing  may  be  assumed  stationary  during  short 
intervals,  either  T =30  msec  for  speech  or  T =100  msec  for  singing. 
The  noise  is  assumed  stationary  during  a longer  period  of  interest. 
During  this  period  of  interest,  therefore,  the  degraded  speech  may  be 
thought  of  as  a concatenation  of  stationary  processes . short  speech 
intervals  degraded  with  additive  stationary  noise . The  noise 
spectrum  generally  has  a common  range  of  frequencies  with  the  speech 
short  interval  spectrum.  To  suppress  the  noise,  therefore,  a Wiener 

filter  can  be  used,  Figure  4-1.  The  Wiener  filter  smoothed  spectrum 
2 

| W(kT ,a))|  is  estimated  for  each  short  interval  T.  The  spectrum  of  a 

2 

rational  filter  H(kT,w)  is  then  matched  to  1/| W(kT,w) | , using 
Autocorrelation  Prediction  described  in  sections  3.4  and  3.5. 

Finally,  the  rational  filter  1/H(kT,w)  is  used  to  filter  the 
corresponding  interval  of  the  degraded  speech,  Figure  4-2. 
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n ( t ) 


Figure  4-1.  Short  interval  Wiener  filtering. 

4.2  Estimation  of  Wiener  Filter 
Smoothed  Spectrum 

2 

The  Wiener  filter  smoothed  spectrum  |w(kT,w)|  for  each  interval 
T is  estimated  from  the  following  equation  * 

W(kT,w)  = 1 - -s (4-1) 

$S+N(kT,U)^ 

A 

where  <J>  (00)  is  an  estimate  for  the  smoothed  spectrum  of  the 
N 

A 

stationary  noise  and  ( kT ,0) ) is  an  estimate  for  the  smoothed 

spectrum  of  the  degraded  speech  during  the  k-th  interval. 

Since  we  are  interested  in  filtering  the  gross  structure  of 
degraded  speech  and  leaving  alone  the  fine  structure,  mainly  caused 
by  the  speech  excitation  function,  the  finite  Linear  Prediction  (LP) 
spectrum  is  used  rather  than  the  standard  Fourier  Spectrum.  Thus, 
<I,j.+N(kT,o')  is  obtained  by  averaging  the  short-time  LP  spectrum  of  the 
degraded  speech  during  the  k-th  interval.  To  estimate  the  noise 


* The  equation  for  Wiener  filter  spectrum  is  derived  in  section  5.5. 
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spectrum,  speechless  portions  of  the  degraded  speech  are  used.  Thus, 

<J>  (u>)  is  obtained  by  averaging  the  short-time  LP  spectrum  during 
N 

these  speechless  portions. 

4 . 3 Pole-Zero  Modeling  Q_f  W ie.ner 
Filter  Smoothed  Spectrum 

Equation  (4.1)  shows  that  the  Wiener  filter  smoothed  spectrum 
|w(kT,w)|2  has  flat  pas3  bands  for  those  frequencies  where  the  speech 
dominates  the  noise,  and  has  deep  stop  bands  at  those  frequencies 
where  the  noise  dominates  the  speech,  Figure  4-3.  In  the  case  where 
speech  dominates  the  noise  at  most  of  the  frequencies,  the  above 
remark  suggests  modeling  |W(kT,w)|2  at  the  stop  bands  first.  This 
can  be  accomplished  by  matching  the  spectrum  of  a rational  filter 
H(kT,w)  to  1/ | W(kT,co) | , using  Autocorrelation  Prediction  described 
in  Section  3-4.  Then,  the  rational  filter  1/H(kT,io)  is  used  to 
filter  the  k-th  interval  of  the  degraded  speech,  Figure  4-2.  In  this 
way,  the  parameters  of  the  model  H(kT,w)  are  reduced  in  number 
considerably  compared  with  an  All-Pole  match  to  |w(kT,co)|  . Figure 
4-4  compares  a Pole-Zero  model  and  an  All-Pole  model  spectral  match 
to  the  same  Wiener  filter  spectrum  in  Figure  4-3.  Note  that  with 
equal  number  of  parameters,  the  Pole-Zero  model  matches  the  deep 
valleys  much  more  accurately  than  the  All-Pole  model. 

Since  H(kT,w)  is  minimum  phase,  the  stability  of  the  rational 
filter  1/H(kT,ui)  is  also  guaranteed. 

4 . 4 Implementation  And  Results 

The  rational  filter  1/H(kT,w)  is  identified  by  two  sets  of 
parcor  parameters  obtained  from  autocorrelation  prediction.  The 
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FREQUENCY  (HZ) 

Figure  4-3.  Wiener  filter  smoothed  spectrum  for  short  interval. 


(a) 


(b) 


Figure  4-4.  (a)  Wiener  smoothed  spectrum  superimposed  with 

30-Pole  and  30-Zero  model  spectrum  (Autocorrelation 
Prediction) . 

(b)  Wiener  smoothed  spectrum  superimposed  with 
60-Pole  model  spectrum  (Linear  Prediction). 
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1/H(kT,w)  is  realized  as  the  cascade  of  two  lattice  filters  described 
in  Section  3-6.  To  improve  the  filtering  process,  both  sets  of  the 
parcor  parameters  are  interpolated  for  each  sample  of  the  input  to 
the  filter. 

In  theory,  the  value  of  W(kT,w)  given  by  (4.1)  is  non-negative 
for  all  frequencies.  In  practice,  however,  W(kT,w)  occasionally 
becomes  negative  at  some  frequencies.  This  is  so  because  of  the 
following:  ^^(kT.w)  is  an  estimate  for  the  speech  smoothed  spectrum 
plus  the  noise  smoothed  spectrum.  On  the  other  hand,  the  interval  T 

(30  msec.)  during  which  the  speech  is  stationary  is  not  long  enough 

* 

to  give  a close  estimate  of  the  noise  smoothed  spectrum.  This  makes 

A 

the  4's+jij( kT , w)  a rough  estimate  and  consequently  W(kT,w),  given  by 
equation  (4-1),  becomes  negative  at  some  frequencies.  To’ ’correct  for 
the  above  problem,  using  the  function  depicted  in  Figure  4-5,  W(kT,w) 
is  tailored  to  positive  values  at  those  frequencies  where  it  is 
negative . 

The  above  Dynamic  Pole-Zero  Filtering  process  was  applied  to 
speech  degraded  with  stationary  additive  colored  noise.  Spectrograms 
of  corresponding  portions  of  the  degraded  speech,  filtered  speech  and 
the  original  clean  speech  are  shown  in  Figure  4-6  for  comparison. 
The  process  was  also  used  to  suppress  helicopter  noise  in  a recording 
of  pilot  s speech  degraded  by  helicopter  noise.  Spectrograms  of 
corresponding  portions  of  degraded  pilot's  speech  and  the  filtered 
speech  are  3hown  in  Figure  4-7.  Figure  4-7a  shows  the  spectrogram  of 
the  helicopter  noise  from  a portion  of  the  recording  where  the  pilot 


was  silent. 
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COMPUTED  MAGNITUDE 


Figure  4-5.  Tailoring  function  (after  T.  Petersen). 

The  Linear  Predition  (LP)  spectrum  of  the  helicopter  noise  is  also 
shown  in  Figure  4-8.  Note  the  harmonic  structure  in  the  helicopter 
noise  depicted  in  both  Figure  4-7a  and  Figure  4-8.  Finally, 
corresponding  portions  of  time  domain  pilot's  speech  with  background 
helicopter  noise  and  the  filtered  speech  are  shown  in  Figure  4-9. 


(b) 


(c) 


Figure  4- 


Corresponding  spectrograms  of  : 

(a)  Input  degraded  speech;  clean  speech  plus  stationary 
colored  noise. 

(b)  Filtered  speech;  output  of  the  Dynamic  Pole-Zero 
Filtering  process. 

(c)  Clean  speech. 

All  three  spectrograms  have  been  scaled  6 dB/oct 
above  400  Hz. 


(b 


(c 


Figure  4-7 


. Spectrograms  of  : 

(a ) Hel icopter  noise. 

(b)  Pilot's  speech  degraded  with  background  helicopter  noise. 

(c)  Filtered  pilot's  speech;  output  of  the  Dynamic  Pole-Zero 
Fi 1 terinq  process . 

All  three  spectrograms  have  been  scaled  6 dB/oct  above 
400  Hz. 
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FREQUENCY  (HZ) 

Figure  4-8.  Linear  Prediction  (LP)  spectrum  of  the  helicopter  noise. 
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Figure  4-9.  Corresponding  portions  of  : 

(a)  Input  pilot's  speech  degraded  with  background 
helicopter  noise. 

(b)  Filtered  pilot's  speech;  output  of  the  Dynamic  Pole-Zero 
Filtering  process. 


CHAPTER  5 


POLE-ZERO  VOCODER  (PZV) 


5 . 1 Introduction 

The  digital  speech  production  model  capitalizing  on  apriori 
information  on  the  structure  of  the  speech  mechanism  and  the  speech 
waveform  is  described.  The  speech  production  model  is  further 
approximated  with  a limited  number  of  parameters,  including  zero  as 
well  as  pole  parameters.  Using  this  parametric  representation  of  the 
speech  production  model,  a so-called  Pole-Zero  Vocoder  (PZV)  is 
devised  for  analysis  and  synthesis  of  clean  speech.  To  code  speech 
degraded  with  stationary  additive  colored  noise,  the  PZV  is  further 
modified  to  account  for  the  noise.  The  PZV,  simulated  on  a computer, 
was  used  successfully  in  analysis  and  synthesis  of  clean  speech 
Finally,  some  pilot  experiments  revealing  the  potential  of  the 
modified  PZV  in  coding  speech  degraded  with  stationary  additive 
colored  noise  is  presented. 

5 . 2 Speech  Production  Model 

For  voiced  sounds,  the  time  varying  filter  model  in  Figure  5-1 
represents  the  effect  of  the  glottis  waveform,  the  vocal  tract,  the 
acoustic  coupling  of  the  nasal  tract  and  the  radiation.  The 
excitation  function  is  a train  of  unit  samples  with  the  same 
frequency  as  the  pitch.  In  contrast,  for  unvoiced  sounds  the  time 
varying  filter  model  represents  the  effect  of  the  vocal  tract  and  the 
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radiation,  and  the  excitation  function  becomes  zero  mean  white  noise 
2 

with  the  variance  a =1. 


Figure  5-1.  Digital  speech  production  model. 

The  time-varying  filter  model  is  considered  to  be  a time-invariant 
filter  during  short  periods  of  time  (10-30  msec.).  This  is  a 
reasonable  assumption,  because  the  dynamics  of  the  articulatary 
configuration  are  slow  due  to  the  inertia  of  muscle  controlled  jaw, 
tongue,  and  lip  movements.  The  finite-time  glottis  waveform,  the 
acoustic  coupling  of  the  nasal  tract  and  the  radiation  cause  the 
transfer  function  of  the  filter  model  to  have  zeros  as  well  the  usual 
poles , [ 1 4 ] . 

The  above  arguments  imply  that  representing  the  filter  model 
transfer  function  by  a finite  dimension  pole-zero  model  [6]  rather 
than  a finite  dimension  all-pole  model  [5]  may  improve  the  quality  of 
the  synthesized  speech  and/or  reduce  the  representation  parameters. 
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5.3  Autocorrelation  Function  of  the  Filter  Model 

We  show  how  the  short-time  autocorrelation  function  of  the 
speech  gives  an  estimate  for  the  autocorrelation  function  of  the 
filter  model  in  Figure  5-1.  The  short-time  autocorrelation  function 
of  the  speech  is  defined  as  the  autocorrelation  function  of  windowed 
speech  using  a smoothed  window  of  proper  length  (20-30  msec.) 

Recalling  the  speech  production  model  in  Figure  5-1,  we  denote 
the  excitation  function,  the  train  of  pulses  or  the  white  noise,  by 
p(n)  and  the  scaled  impulse  response  of  the  filter  model  by  v(n). 
The  speech  s(n)  is,  therefore,  expressed  as: 

s (n)  = p (n)  © v (n)  . (5-1) 


To  obtain  a stationary  segment  of  speech,  the  speech  s(n)  is  weighted 
by  some  smoothed  window  w(n)  of  proper  length  (20-30  msec.),  i.e. 


s 

w 


(n) 


A 


s (n) 


w(n) 


[p(n)  ® v(n)  ] w(n) 


(5-2) 


Assuming  the  window  w(n)  is  smooth  during  the  effective  duration  of 
the  impulse  response  v(n) , the  equation  (5-2)  can  be  approximated  as: 


s (n)  - [p (n)  w(n) ] ® v(n) , (5-3a) 

w 

or 

(n)  ~ Pw(n)©v(n),  (5-3b) 

where  Pw(n)  I3  the  weighted  excitation  function  defined  as: 
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p (n)  = p(n)  w(n).  (5-4) 

w 

Short-time  Autocorrelation  Function  of  Speech  Using  (5-3b),  the 

speech  short-time  spectrum  S (z )S  (i ) is  approximated  by: 

w w ^ 

S (z)S  (1/ z)  = V(z)V(l/z)  P (z)P  (1/z).  (5-5) 

w w w w 

The  equation  (5-5)  shows  that  the  speech  short-time  spectrum 
approximates  the  filter  model  spectrum  multiplied  by  the  spectrum  of 
the  windowed  excitation  function.  From  (5-5),  the  speech  short-time 
autocorrelation  function  R(k)  is  approximately  equal  to 

R(k)  = R (k)  © R (k),  (5-6) 

v rw 

where  Ry(k)  and  Rp  (k)  are  autocorrelation  functions  of  the  filter 
w 

model  and  the  windowed  excitation  function,  respectively. 

Now  we  find  the  Rp  (k)  for  both  voiced  and  unvoiced  speech.  In 
w 

the  case  of  voiced  speech,  the  excitation  function  p ( n ) is  a train  of 
pulses  with  the  same  period  T as  the  pitch  period.  Thus  using  (5-4), 

the  weighted  excitation  function  p (n)  takes  the  form  of  a weighted 

w 

train  of  pulses,  i.e. 


+oo 

Pw(n)  = w(n)  j 6(n-(*T) 

f =s— OO 


n = 0,  +T , . . . 
otherwise . 


(5-7) 


The  autocorrelation  function  of  Pw(n),  therefore,  becomes 


Thus  for  voiced  speech,  the  short-time  autocorrelation  function  R(k) 

is  approximately  equal  to  the  filter  model  autocorrelation  function 

Ry(k)  convolved  with  the  symmetric  decaying  train  of  pulses  R (k), 

Pw 

having  the  same  period  T as  the  pitch  period,  see  Figure  5-2a. 
Hence,  the  R ( k ) is  an  estimate  of  Rv(k)  for  short  time-lags.  For 
unvoiced  speech,  the  excitation  function  p(n)  is  white  noise. 
Assuming  the  smoothed  window  w(n)  is  long  enough,  then  the 
autocorrelation  function  of  Pw(n)  approximates  the  autocorrelation 
function  of  white  noise  namely  ,<5(n) . Thus  we  have 


R (k)  = 6 (k ). 
P 
w 


(5-9) 


As  a result,  using  (5-6)  the  short-time  autocorrelation  function  for 
unvoiced  speech  is  an  estimate  for  the  filter  model  autocorrelation 
function  Rv(k),  see  Figure  5-2b. 


Each  short  segment  of  speech  can  be  represented  by  a 
Voiced/Unvoiced  decision,  the  pitch  period,  if  the  decision  is 
voiced,  and  the  filter  model  transfer  function.  To  be  able  to 
represent  the  short  segments  of  speech  with  a limited  number  of 


parameters,  the  filter  model  in  Figure  5-1  is  represented  by  a 


TIME  LAG 

(a) 


HMC't 

• ‘ * i Mica 


TIME  LAG 

(b) 

Figure  5-2.  (a)  Short-time  autocorrelation  function  of 

voiced  sound  [0]. 

(b)  Short-time  autocorrelation  function  of 
unvoiced  sound  [P]. 
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parametric  model.  Recalling  the  discussion  in  section  5.2,  the 
parametric  model  is  chosen  to  be  the  Pole-Zero  model  H(z)  given  by 
equation  (3-3).  Furthermore,  since  the  numbers  of  the  spectral  poles 
and  the  spectral  zeros  in  the  filter  model  spectrum  are  different  for 
different  segments  of  speech,  then  full  advantage  of  the  Pole-Zero 
modeling  is  gained  by  allowing  the  orders  M and  L to  be  dynamic. 

The  overall  block  diagram  of  the  so-called  Pole-Zero 
Vocoder  (PZV)  is  shown  in  Figure  5-3-  The  Voiced/llnvoiced  decision 
and  the  pitch  period  are  extracted  using  any  of  a variety  of  existing 
methods  [28].  The  parameters  of  the  Pole-Zero  model  H(z)  are 
obtained  by  applying  Autocorrelation  Prediction,  Figure  3-2,  to  the 
short-time  autocorrelation  function  of  the  speech.  Since  the  speech 
short-time  autocorrelation  function  is  an  estimate  for  the  filter 
model  autocorrelation  function,  then  the  spectrum  of  the  resulting 
Pole-Zero  model  approximates  that  of  the  filter  model.  In  applying 
Autocorrelation  Prediction  to  the  speech  short-time  autocorrelation 
function,  the  parsimoneous  (most  economical)  orders  Mp  and  Lp  for 
each  segment  may  be  obtained  using  the  methods  described  in  Section 
2-5.  Finally,  to  synthesize  speech  using  the  PZV  depicted  in  Figure 
5-3,  we  update  the  Voiced/Unvoiced  decision,  the  pitch  period  T and 
the  Pole-Zero  model  parameters  for  each  segment  of  speech. 

5 • 5 Pole-Zero  Analysis-Synthesis  of  Speech 

In  Presence  of  Additive  Stationary  Noise 

The  idea  is  to  find  the  optimal  estimate  of  the  speech  spectrum 
for  each  short  interval.  The  autocorrelation  function  of  this 


optimal  estimate  is  then  matched  with  that  of  a Pole-Zero  model  using 


Figure  5-3.  Pole-Zero  Vocoder  (PZV)  overall  block  diagram. 
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Autocorrelation  Prediction  described  in  section  3-4.  We  make  the 
assTimption  that  the  background  noise  is  stationary  during  "a  long 
period  of  interest,  while  the  speech  is  stationary  only  in  3hort 
intervals  (30  msec). 

The  optimum  linear  filter,  namely  the  Wiener  filter,  is  designed 
for  each  short  interval  of  noisy  speech,  Figure  5-4. 


n(t) 


Figure  5-4.  Wiener  filtering 
The  transfer  function  of  the  Wiener  filter  is 


H (60  ) 


'sx(a)) 
*xx(u)  ■ 


(5-10) 


Assuming  the  noise  and  speech  are  uncorrelated,  then  (5-10)  takes  the 
form 

A 


H(w) 


*ss(<»>) 


*SS(U))  + Wu) 


1 + 


W“>’ 

♦ss0-) 


(5-11) 


or 
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, . ' '*xx(u)  - $nn(w)  VU) 

H (oj)  = : ; ; “ l 


V“> 


V“> ' 


(5-12) 


The  noise  spectrum  ^^(w)  can  be  estimated  using  speechless  portions 
from  the  noisy  speech  record.  From  Figure  5-4,  the  optimal  estimate 
of  the  short  interval  spectrum  is 


= $>  v(u>)|H(u))|2.  (5-13) 

ss  xx 


L Using  (5-11)  or  (5-12)  in  (5-13)  gives 


(<»)) 

SS 


4> 


1 + 


ss (co) 
*ss(w) 


(5-14)  . 


or 


Vw)  = *xx(w) 


1 - 


*xx(w) 


(5-15) 


Equation  (5-14)  shows  how  the  ratio  of  the  optimal  estimate  $^(w)  to 

SS 

the  speech  spectrum  $ss(w)  is  related  to  the  power  ratio 


$.,»,(<*>)/$„_ (w)  of  the  background  noise  to  the  speech  at  different 
NN  j SS 
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frequencies.  The  optimal  estimate,  however,  is  computed  from  (5-15), 

where  <J>  (w)  is  the  noisy  speech  spectrum  for  the  short  interval  (30 

XX 

msec)  and  <J>  (to)  is  the  noise  spectrum  for  the  period  of  interest. 

NN 

An  estimate  for  is  obtained  by  averaging  spectra  or  LP 

spectra  for  short  overlapping  segments  of  noisy  speech.  An  estimate 

for  ^NN(W)  is  computed  by  averaging  the  short-time  spectra  or  LP 

spectra  during  speechless  segments  of  noisy  speech.  Finally,  the 

Inverse  Fourier  Transform  of  the  optimal  estimate  <Jv~((jj)  gives  the 

SS 

autocorrelation  function  (k),  i.e. 


R(k) 


<J>^(<jd)e  -1kwdu). 
SS 


(5-16) 


The  overall  block  diagram  for  the  analysis  and  synthesis  of 
speech  in  the  presence  of  additive  stationary  noise  is  depicted  in 
Figure  5-5.  Similar  to  clean  speech  analysis,  applying 

Autocorrelation  Prediction  to  the  R(k)  gives  the  parameters  of  the 
Pole-Zero  model,  the  autocorrelation  function  of  which  approximates 
the  R(k).  The  Voiced/  Unvoiced  decision  and  the  pitch  period 
extraction  for  the  short  segments  of  speech,  becomes  more  complicated 
in  the  presence  of  noise  [ ] . The  synthesis  remains  the  same  as  the 
one  in  PZV. 

5 . 6 Implementation  and  Results 

The  Pole-Zero  Vocoder  (PZV)  based  on  Autocorrelation  Prediction 
was  simulated  on  a PDP-10  computer.  The  parameters  of  a fixed  order 
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Pole-Zero  model  are  estimated  by  applying  Autocorrelation  Prediction 
to  the  speech  short-time  autocorrelation  function.  The  Pole-Zero 
model  is  realized  by  a lattice  filter  augmented  with  tap  parameters 
[17].  The  gain  G for  the  lattice  filter  is  computed  such  that  the 
corresponding  synthesized  speech  and  the  natural  speech  have  equal 
energy  [5,8]. 

The  Voiced/Unvoiced  decision  and  the  pitch  period  is  extracted 
using  an  implementation  of  the  SIFT  algorithm  [28]  by  Boll  [9].  The 
analysis-synthesis  is  quasi-pitch  synchronous  [9].  The  augmented 
lattice  filter  parameters  are  linearly  interpolated  between  two 
successive  analysis  segments. 

This  vocoder  was  used  for  analysis  and  synthesis  of  passages  of 
natural  speech.  For  proper  selection  of  the  fixed  orders  L and  M, 
informal  hearing  tests  show  some  improvement  of  the  synthesized 
speech  generated  by  the  Pole-Zero  Vocoder  over  that  generated  by  the 
All-Pole  vocoder  having  comparable  number  of  parameters.  The 
improvement  is  more  noticeable  when  the  corresponding  natural  speech 
has  more  nasalized  sounds.  The  Pole-Zero  synthesized  speech  has  less 
"ringing"  quality  than  the  corresponding  All-Pole  synthesized  speech. 
Figure  >-6  shows  corresponding  segments  of  natural  Pole-Zero 
synthesized  and  All-Pole  synthesized  speech. 

A pilot  experiment  to  test  the  performance  of  the  Pole-Zero 
analysis-synthesis  of  speech  in  presence  of  additive  colored  noise 
was  implemented  as  follows:  Using  the  dynamic  pole-zero  filtering 
process  described  in  Chapeter  4,  the  noisy  speech  was  first  filtered 
and  then  the  Pole-Zero  Vocoder  (PZV)  was  applied  to  the  resulting 
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filtered  speech.  Figure  5-7  shows  the  spectrograms  of  corresponding 
portions  of  vocoded  noisy  speech,  vocoded  filtered  speech  and  vocoded 
clean  speech.  .Informal  hearing  tests  reveal  an  improvement  in  the 
vocoded  filtered  speech  over  the  vocoded  noisy  speech.  The  Pole-Zero 
coding  of  speech  in  the  presence  of  noise,  described  in  section  5.5, 
which  combines  the  Pole-Zero  filtering  of  noisy  speech  and  the 
Pole-Zero  coding  of  clean  speech,  seems  to  have  the  potential  to 
generate  improved  coded  speech. 


AMPLITUDE 


Figure  5-5.  Corresponding  segments  of  : 

(a)  Voiced  natural  SDeech. 

(b)  Pole-Zero  synthesized  speech  (12P  + 4Z). 

(c)  All-Pole  synthesized  speech  (16P). 


(a ) 


(b) 


(c) 


Figure  5-7.  Spectrograms  of  corresponding  portions  of  : 

(a)  Vocoded  filtered  speech. 

(b)  Vftcoded  clean  speech. 

(c)  Vocoded  noisy  speech. 

All  three  spectrocrams  have  been  scaled  6 dB/oct 
above  400  Hz. 


CONCLUSIONS 


6 . 1 Review 

The  goal  of  this  research  has  been  to  develop  new  Pole-Zero 
modeling  techniques  and  apply  them  to  speech  processing.  This  goal 
has  been  accomplished.  Two  Pole-Zero  modeling  techniques, 
Autocorrelation  Partial  Realization  (APR)  and  Autocorrelation 
Prediction  (AP),  have  been  developed  and  their  theories  were 
established.  APR,  using  partially  linear  operations,  identifies  the 
Pole-Zero  model  whose  short  time-lag  autocorrelations  exactly  match 
those  of  a given  spectrum.  in  contrast  AP,  using  only  linear 
operations,  identifies  the  Pole-Zero  model  whose  short  time-lag 
autocorrelations  closely  approximate  those  of  a given  spectrum. 
Neither  of  them  uses  Fourier  Transformation,  but  fast  recursive 
and/or  iterative  algorithms  to  estimate  the  model  parameters.  APR 
and  AP  have  been  compared  and  the  properties  of  the  Pole-Zero  models 
identified  by  them  were  discussed.  It  has  been  shown  that  the 
Pole-Zero  model,  identified  by  Autocorrelation  Prediction,  has 
advantages  over  the  All-Pole  model,  identified  by  Linear 
Prediction  (LP),  when  the  the  envelope  of  the  given  spectrum  has  deep 
valleys.  A cascade  of  two  lattice  filters  has  been  proposed  as  a 
realization  of  the  Pole-Zero  models  identified  by  AP  or  APR. 

A dynamic  filtering  process,  based  on  Wiener  filtering  and 
Autocorrelation  Prediction  has  been  developed  and  implemented  to 
suppress  background  noise  from  degraded  speech.  Using  AP,  rather 
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than  LP,  to  model  the  estimated  Wiener  filter  spectrum,  has  improved 
the  performance  of  this  dynamic  filtering  process.  Moreover,  a 
Pole-Zero  Vocoder  (PZV)  based  on  AP  has  been  developed  and 
implemented.  PZV  becomes  an  All-Pole  Vocoder  (Linear  Predictive 
Vocoder)  on  one  extreme  and  an  All-Zero  Vocoder  (different  from  the 
Homomorphic  Vocoder  but  similar  in  quality  of  the  synthesized  speech) 
on  the  other  extreme.  For  proper  selection  of  the  Pole-Zero  model 
orders  (M  is  roughly  three  times  L),  informal  tests  has  shown  that 
PZV  generates  more  "natural  sounding"  synthesized  speech  than  the 
other  two  extremes. 

6 . 2 Future  Research 

To  use  the  Pole-Zero  model  obtained  from  APR  in  applications 
other  than  spectral  matching,  like  filtering,  the  stability  of  the 
model  should  be  resolved.  This  area  needs  futher  investigation. 

To  obtain  a low  bit  rate  high  quality  vocoder,  the  orders  of  the 
Pole-Zero  model  in  PZV  can  be  made  dynamic;  preferably  equal  to  the 
parsimonious  orders  M and  L for  each  frame  of  anlysis. 

Modifying  the  PZV  to  account  explicitly  for  the  background  noise 
of  degraded  speech  seems  promising  when  the  noise  is  stationary  and 
only  the  degraded  speech  is  available.  The  restriction  of  the 
stationarity  of  the  background  noise  may  be  relaxed  if  a correlated 
version  of  the  background  noise,  recorded  by  another  microphone  away 
from  the  speaker,  is  also  available. 

Finally,  the  applications  of  Autocorrelation  Prediction,  as  a 
Pole-Zero  modeling  technique,  in  areas  other  than  speech  processing 
are  open  for  further  research. 


APPENDIX  A 


DURBIN  AND  PARCOR  RECURSIVE  ALGORITHMS 

A system  of  linear  equations  can  be  solved  recursively  by  the 
Bordering  method  [l3].  When  the  coefficient  matrix  of  the  system  of 
linear  equations  is  a Symmetric  Toeplitz  matrix,  then  the  Bordering 
method  is  simplified  to  Levinson's  algorithm  [42].  If  the  constant 
vector  in  the  right  hand  side  of  the  system  of  linear  equations  is 
also  of  the  following  form,  equation  Al,  then  Levinson's  algorithm  is 
further  simplified  to  Durbin's  algorithm  [12,  27]. 
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Durbin's  Recursive  Algorithm 

The  system  of  linear  equations  (Al ) can  be  solved  recursively  by 


Durbin's  algorithm: 
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1) 

E0  = R0’  j = °’ 

2) 

j = j + 1, 

3) 

K.  = - R.  + l afj 

J Lj  iti  1 

V<]  /v- 

4) 

a«>  = K., 
J J 

5) 

a9^  = + K. 

J i J 

•ft”. 

6) 

Ej  - a - k*)^. 

7) 

If  j = M stop;  otherwise  go  to  2 

When  the  algorithm  stops,  the  solution  to  the  system  of  equation 
(Al)  is: 

(M) 

a.  = a^“' , 1 < i < M . (A3) 


The  corresponding  by-product  parameters  K^,  for  1 £ i £ M are 
referred  to  as  the  reflection  coefficients,  partial  correlations  or 
parcor  parameters. 

If  the  function  represents  an  autocorrelation  function,  then 
the  coefficient  matrix  in  (Al)  is  positive  definite  [34],  In 
this  case,  it  can  be  shown  that  the  parcor  parameters  have  the 
property 


< 1, 


1 < i < M, 


(A4 ) 
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M 

and  the  polynomial  A(z)  = 1 + £ a.  z is  minimum  phase.  That  is, 

i=l  1 

all  the  roots  of  A(z)  are  inside  the  unit  circle  [ 18 ] . 

There  is  a unique  set  of  parcor  parameters  for  the  set  of 

predictor  coefficients  (a.J  and  vice  versa.  To  compute  directly  the 
corresponding  parcor  parameters  {k^lfrom  the  predictor  coefficients 
{a^},  a recursive  algorithm  is  derived.  Substituting  j-1  for  i in 
step  5 of  Durbin's  algorithm  and  rearranging  the  terms  gives 


a 


(j-1) 

J-i 


K.  a 
1 


(j-1) 


1 < i < j-1. 


(A5) 


Again  substituting  for  in  step  5 

from  (A5)  and  some  simplifications  results  in 


of  Durbin's  algorithm 


afj-1) 

1 


K. 

1 


(1-Kj ) , 


1 1 i < j-1. 


(A6) 


Thus  from  step  A in  Durbin's  algorithm,  (A3)  and  (A6)  we  obtain  the 
sought  for  parcor  recursive  algorithm. 


1)  J = M,  aJM)  = a.. 


1 < i < M 


2) 


K = a(j) 

Kj  j ’ 


3) 

A) 

5) 


If  j = 1 stop 


,<>»  • [.«>  - Kj  .0>J  /„- «*,. 


1 < i < j-1  , 


J - j - 1 , 


6) 


Go  to  2 . 


(A7 ) 


APPENDIX  B 


PADE  APPROXIMATION 


A function  represented  by  a one-sided  power  series  is  approxi- 
mated with  a rational  function  using  the  Pade  approximation. 

Let  X(z)  be  a function  represented  by  the  following  one-sided 
power  series 


X(z)  = l x z k,  for  | z | > r . 

k=0 


Also,  consider  the  rational  function 

-i 


L 

l Ci- 

y/z\  = C(z)  = *-~fl 

X(  } D(z)  M 

y d . z 

u 1 
i=0 


d0  1’ 

dM*  °’ 
L < M, 


(Bl) 


(B2-a) 


whose  power  series  representation  is 

OO 

X(z)  = l £ z for  | z | > r„  . (B2— b ) 

k=0 

4 

From  where  (B2-a)  and  (B2-b) , the  coefficient  x^  is  equal  to 
L M 

xk  * l ~ l dl  \-i  k > 0 . (B3) 

K i-0  i-1 


( 
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To  approximate  the  X(z)  by  X(z) , the 

N = L + M + 1,  (RM 

unknown  coefficients  {d^}  and  {c^}  are  computed  such  that  the  first  N 

A 

terms  of  the  x(z)  and  X(z)  power  series  representations  are  equal, 
i . e. 
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= I 


i=0 


d , x,  . 
1 k-i 


d = 1 , 

o ’ 

0 < k < L 


(B7-b) 


M 

where  the  matrix  (xT  , , . , (B7-a)  is  non-symmetr ic  toeplitz.  The 

L+i-j ) i 

Trench  recursive  algorithm  [39,  Appendix  C]  can  be  used  solve  the 
system  of  linear  equations  (B7-a).  After  computing  the  {d^}  the 
coefficients  {c^}  are  obtained  simply  by  performing  the  summation  in 
(B7-b) . 


Pade  Approximation  Error 


Uaing  (Bl),  (B2-b)  and  (B5) , the  Pade  approximation  error.  Figure 
B-^.,  is  defined  by  the  following  power  series 


Figure  B-l.  Pade  approximation  error 


81 


E(z)  = I e.  z k = x(z)  - x(z)  = l (x  - x )z  k 
k=0  k=0  K 


= [ (xk  - xk)z 


-k 


k=N 


|z|  > Ma x ( r j , rp 


(B8-a) 


A 

Using  (B5)  and  substituting  for  x^  from  (B3)  in  (B8-a)  gives 
another  expression  for  E(z),  i.e. 


E (z ) 


CO 


oo 


l 

k=N 


i=0 


d. 

1 


(B8-b) 


Relation  (B8-a)  shows  that,  the  first  discrepancy  between  the 
corresponding  power  series  coefficients  of  the  X(z)  and  X(z)  occurs 
at  k * N.  The  relation  (B8-b)  shows  that  different  errors  E(z)  are 
obtained,  depending  on  the  selection  of  the  orders  L and  M for  a given 
N.  All  of  these  different  errors  E(z),  nevertheless,  have  the  following 
property: 

ek  = 0 , 0 < k £ N-l  . (B9) 

Special  Case 

>.  C (z  ) 

The  Pade  approximation  for  the  special  case  L=M,  and  conse-  ^ • 

quently  N=2M+1,  is  considered  here.  Also,  the  value  of  the  first 
discrepancy  e^  is  expressed  in  terms  of  only  the  given  power  series 

k' 


coefficients  x 
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The  matrix  equation  (B7-a)  for  L=M  takes  the  form 


*M 

XM-1 

X2 

X1 

dl 

XM+1 

XM+1 

*M 

X3 

x2 

d2 

XM+2 

• 

# * 

• • 

* 

— _ 

* 

X2M-2 

X2M-3 

' ’Si 

*m-i 

dM-l 

X2M-1 

X2M-1 

X2M-2 

‘ XMfl 

*M 

dM 

X2M 

(BIO) 


The  non-symraetr ic  toeplitz  matrix  ^)M  in  (BIO)  is  symmetric 


around  its  second  diagonal.  Therefore,  by  rearranging  the  unknown 
coefficients  {dD  of  (BIO)  in  reverse  order,  the  matrix  equation  (BIO) 
takes  the  following  form 


X1  X2  * • Vl  XM 

dM 

’Si+i 

X2  X3  * • "M  Vi 

dM-l 

XM+2 

XM-1  ’Si  ■ • X2M-3  X2M-2 

d2 

X2M-1 

’Si  XM+1  ' ' X2M-2  X2M-1 

dl 

X2M  " 

_ _ 

(Bll) 


M 

The  symmetric  matrix  (x  _^)^  in  (Bll)  is  a Hankel  matrix. 
The  Berlekamp-Massey  algorithm  (Appendix  C)  can  be  used  to  solve  the 
system  of  linear  equations  (Bll).  From  (B8-b)  and  (B5),  the  first 
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discrepancy  e^,  for  the  case  L=M  is  given  by 


eN  C1M+1  di  X2M+l-i 

1=0 


To  compute  e^  directly  in  terms  of  only  the  given  power  series  coef- 
ficients , the  system  of  equations  (Bll)  is  augmented  by  the 

equation  (B12).  e^  is  considered  as  an  unknown.  This  leads  to 


i+j-1 


i+j-1 


where 


lx.,.  , is  the  determinant  of  the  Hankel  matrix  (x.,.  ,)i 

I 1+J-:  1 , ,M+1  1+3  “L  1 


given  in  (Bll)  and  the  x 


Hankel  matrix: 


i+j-i i 


is  the  determinant  of  the  following 


X1  X2 


iVl 


X2  X3  *M+1  XM+2 


*M  XM+1 


X2M-1| X2M 


XM+1  XM+2  ‘ ' X2M  X2M+1 


appendix  c 


\ 


TRENCH  AND  BERLEKAMP-MASSEY  RECURSIVE  ALGORITHMS 


In  using  the  Pade  approximation,  one  encounters  the  system  of 
linear  equations  of  the  form 


— 

■* 

^0  *-1  *-2  • • ^-M+l 

*1 

♦i 

*1  *0  *-1  • * *-M+2 

^2 

^2 

*2  01  *0  • • ^-M+3 

-s- 

— — 

^3 

^M-l  ' M-2  '^M-3  • • '0 

Ip 

(Cl) 


M 

where  the  coefficient  matrix  (0^  . )j  in  (Cl)  is  a non-symmetr ic 
Toeplitz  matrix.  Using  the  Bordering  method  [ 1 3 ] , Trench  [39], 
and  Whittle  [4 1 ] d erive  a fast  recursive  algorithm  to 
solve  simultaneously  (Cl)  and  the  following  system  of  equations 


— — 

— 



^0  h ^2  •••  Vi 

T1! 

♦-1 

♦-1  ^0  *1  •••  *M-2 

n2 

♦-2 

*-2  ^-1  *0  *M-3 

n3 

— 

♦-3 

__^-M+l^_ff4-2 *-M+3  ' V°  _ 

(C  2 ) 
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Trench's  Recursive  Algorithm: 


1) 

> 

o 

II 

o 

J = o. 

2) 

j = j + 1> 

3) 

= - 

r J-i  O-i) 

h + ^ *i  *i-i\  / 

A) 

II 

A-1  • 

5) 

m (j-D  (j)  (j-D 

iKJ'  = iii  + \ b.  n 

i i J j-i 

1 < i < j-1  , 

6) 

= n! 

!j_1)  + nfj)  , 

l J J-i 

1 < i < j-1  , 

7) 

AJ  = (i  - 

- ^ , 

8) 

If  j = M 

stop;  otherwise  go  to  2. 

The  system  of  equations  (Cl)  becomes  identical  to  the  system  of 
equations  (B7-a)  by  replacing  by 

-L  < k < M-l  , 

Otherwise  , (C3) 

and  by  d^  for  1 f.  k <_  M.  Thus  Trench's  recursive  algorithms  can  be 
used  to  solve  the  system  of  equations  (B7-a)  encountered  in  the  Pade 
approximation.  Note  also,  the  system  of  equations  (Cl)  becomes  identi- 
cal to  (C2),  except  for  the  unknown  names,  when  the  coefficient  matrix 
M 

is  symmetric,  i.e.,  = <f)_k 


In  this  case,  the  Trench 
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algorithm  reduces  to  Durbin's  algoritlun. 


1)  1 ■>  D(z)  1 - P(z)  1 - 8 

0 - M 1 - p 0 -►  k . 

2)  If  k = N,  stop;  otherwise  compute 

M 

e = x.  + l d . x,  . . 
k ,=1  i k-i 

3)  If  e = 0,  then  8 + 1 - 8 and  go  to  6). 

A)  If  e / 0 and  2M  > k,  then 

— 1 ~ K 

D(z)  - ep  z P(z)  -*•  D(z) 

8+1-6 
and  go  to  6) . 

5)  If  e / 0 and  2M  <_  k,  then 

D(z)  - T(z)  (temporary  storage  of  D(z)) 

— i — R 

D(z)  - ep  z P(z)  - D(z) 
k + 1 - M - M 
T(z)  - P(z) 
e - p 

1 - 8- 


6) 


k + 1 ♦ k and  return  to  2) . 
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where 

D(z):  current  denominator 

M:  order  of  current  denominator 

p(z):  last  denominator  of  lower  order 

e:  next  discrepancy 

p:  the  discrepancy  corresponding  to  the  last  denomination 

of  lower  order 

k:  current  number  of  matched  terms 

8:  difference  between  current  number  of  matched  terms  and 

that  corresponding  to  last  denominator  of  lower  order. 


The  Berlekamp-Massey  algorithm  [7,  30]  recursively  finds  the 


weights  (d . } of  the  shortest  Linear  Feedback  Shift  Register  (LFSR) 
that  ran  generate  exactly  the  given  sequence  fx.  Figure  C-l. 

K.  1 — U 

The  LFSR  given  in  Figure  C-l  represents  the  rational  function 
M-l 

c ' . 

l 


C’  (z)  = _i=J> 

D(z) 


I c ' • z_  1 


Lv'1-  do = l- 


M 

y d . z 

i=0 


-1  i=0 


(CM 


The  coefficients 


sequence  lx. 


,N-1 


k 0 

the  weights  (d ^ } 


{ c ^ } are  obtained  from  the  first  M values  ! r 
stored  in  the  shift  register  as  initial  va 
using  the  following  relation: 


c 


M 


l 

k=0 


dkXi-k 


0 < i < M-l  . 


Hence  one  can  use  the  Ber lekamp-Masse  i 
nator  coefficients  of  the  least  dimt  -i 
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(C4)  for  which  the  first  N terms  of  its  power  series  representations 
are  given. 

We  are  interested,  however,  in  the  rational  function  of  the 


form 


cl?.) 

D(z) 


M 

l 

i=0 


M 1 

l V 


M-l 


M 1 

l at. 1 

i-0 


c0  + 2 


l C'i 

-1  i-0 


•z'1 


n i 

l diz 

i-0 


= x° + ii  v'1,  do*1. 


(C5) 


rather  than  the  form  (C4).  From  (C5)  and  (C4)  one  concludes  that  the 

denominator  coefficients  {d^} in  (C5)  are  the  feedback  weights  of  the 

N-l 

shortest  LFSR  which  generates  the  sequence  (x  } rather  than  the 
N-l 

sequence  {x^}.  q . 

Therefore,  applying  the  Berlekamp-Massey  algorithm  to  the  sequence 
N— 1 

{ xi } i_ i gives  the  desired  denominator  coefficients  {d ^ . Using  (C5) , 
the  numerator  coefficients  are  obtained  from  the  following  formula 


i 

ci  = I dk  Xl_k  , 0 < i < M . (C6) 

k-0 

Figure  C-2  represents  the  rational  function  defined  in  (C5). 

If  the  given  sequence  was  indeed  generated  by  some 

unknown  rational  function  of  the  form  (C5)  and  the  condition  N 1 2M+1 
holds,  then  the  Berlekamp-Massey  algorithm  detects  the  order  M and 
gives  the  corresponding  denominator  coefficients  {d^}^^.  When  the 
sequence  corresponds  to  some  real  data,  however,  the  Berlekamp- 

Massey  algorithm  applied  to  the  sequence  {x  , recursively  computes 

the  solution  to  the  following  system  of  equations,  (C7),  of  order 
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M = [n/2].  The  [N/2]  means  the  largest  Integer  less  than  or  equal 
to  the  real  number  N/2. 


X1  x2  x3  • • • *M 

1 

1 ^ 

Vi 

X2  x3  X4  * ' ' XHfl 

Vi 

XW-2 

X3  X4  X5  * • ‘ XNH-2 

• • • • 

• • • • 

dM-2 

— “ 

Sl+3 

• • • • 

’Si  Sft-l  *Mf2  • • • X2M-1 

dl 

X2M 

a,  _ 

(C  7 ) 


APPENDIX  D 


NEWTON-RAPHSON  ITERATIVE  ALGORITHM 


To  decompose  the  symmetric  polynomial 


M i 

PU) ' 1 ..  pui  2 • 


into  the  product  of  the  form 


/ M i „ , M , . M , M— I i | . 

KmWM  ■ Uo  bi*  ) Li  b,z ) ■ tL„(  i biVu)- 


the  polynomial  P(z)  should  be  non-negative  on  the  unit  circle  in  z- 
plane.  Assuming  P(z)  has  this  property,  an  iterative  algorithm  based 
on  Newton-Raphson  method  is  formulated  to  compute  the  polynomial  B(z) 
such  that 


P(z)  - 


) - J_M(p,ir|_o  Vj+ui)*'1  -0- 


Consider  the  vector  function  (f^)  whose  i-th  element  is  defined 


[i ' pi  - bjbj«  • 0 i 1 i " 


<?2 

Computing  the  B(z)  such  that  the  relation  (D3)  holds,  is  equiva- 
lent to  finding  the  vector 


(D5) 


for  which  the  vector  function  (f^  has  zero  value.  Therefore,  the 
Newton-Raphson  iterative  method  can  be  used  to  find  the  vector  (b^) 
and  as  a result,  the  desired  polynomial  B(z).  The  iterative  Newton- 
Raphson  formula  for  the  vector  function  (f^)  takes  the  form 


3(f,) 


<«t>  * t<bi>t+I  - <>>1)']  5(5^  ‘ 0 


(D6) 


where  (b^)*"  is  the  estimate  for  the  desired  (b^)  at  the  t-th 
iteration. 

A a(fi> 

From  (D2) , the  matrix  T ■ ^-/u  \ is  equal  to 


3(fi>  l*Jt\ 

30^)  “ l 3b J J 


3(b  ) 

“ 1 

b0  bl  ' * ’ bM 


bl  b2 


M 


• ■ ' o 


M 


b0  bl 


O 


M 

bM-l| 


(D7) 


Using  (D7)  in  (D6)  and  rearranging  the  terms,  results  in 


(b1)t4’1  - (b1)t  - t"1  (f4). 


(D8) 


Q3 


The  relations  (DA)  and  (D7)  and  (D8)  provide  the  iterative 
algorithm  for  computing  the  desired  vector  (b^). 

Assuming  pQ>0,  the  starting  vector  (b^0  can  be  chosen  to  be 


(D9) 

The  iteration  is  considered  to  have  converged  when  for  some  pre- 
scribed value  e the  following  inequalities  hold 

fj  < £ , 0 < i < H . (DIO) 

It  can  be  shown  that  the  convergence  is  of  second  order  and  the 
polynomial  Bt(z)  corresponding  to  the  vector  (bi)t,  having  the  (b^0 
as  starting  vector,  is  minimum  phase  [43]. 
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